Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/MultiSpecies/DGMultiSpecies.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Compressible multi-species flow using discontinuous Galerkin
9 : : finite elements
10 : : \details This file implements calls to the physics operators governing
11 : : compressible multi-species flow using discontinuous Galerkin discretization.
12 : : */
13 : : // *****************************************************************************
14 : : #ifndef DGMultiSpecies_h
15 : : #define DGMultiSpecies_h
16 : :
17 : : #include <cmath>
18 : : #include <algorithm>
19 : : #include <unordered_set>
20 : : #include <map>
21 : : #include <array>
22 : :
23 : : #include "Macro.hpp"
24 : : #include "Exception.hpp"
25 : : #include "Vector.hpp"
26 : : #include "ContainerUtil.hpp"
27 : : #include "UnsMesh.hpp"
28 : : #include "Inciter/InputDeck/InputDeck.hpp"
29 : : #include "Integrate/Basis.hpp"
30 : : #include "Integrate/Quadrature.hpp"
31 : : #include "Integrate/Initialize.hpp"
32 : : #include "Integrate/Mass.hpp"
33 : : #include "Integrate/Surface.hpp"
34 : : #include "Integrate/Boundary.hpp"
35 : : #include "Integrate/Volume.hpp"
36 : : #include "Integrate/Source.hpp"
37 : : #include "Integrate/SolidTerms.hpp"
38 : : #include "RiemannChoice.hpp"
39 : : #include "MultiSpecies/MultiSpeciesIndexing.hpp"
40 : : #include "Reconstruction.hpp"
41 : : #include "Limiter.hpp"
42 : : #include "Problem/FieldOutput.hpp"
43 : : #include "Problem/BoxInitialization.hpp"
44 : : #include "PrefIndicator.hpp"
45 : : #include "MultiSpecies/BCFunctions.hpp"
46 : : #include "MultiSpecies/MiscMultiSpeciesFns.hpp"
47 : : #include "EoS/GetMatProp.hpp"
48 : :
49 : : namespace inciter {
50 : :
51 : : extern ctr::InputDeck g_inputdeck;
52 : :
53 : : namespace dg {
54 : :
55 : : //! \brief MultiSpecies used polymorphically with tk::DGPDE
56 : : //! \details The template arguments specify policies and are used to configure
57 : : //! the behavior of the class. The policies are:
58 : : //! - Physics - physics configuration, see PDE/MultiSpecies/Physics.h
59 : : //! - Problem - problem configuration, see PDE/MultiSpecies/Problem.h
60 : : //! \note The default physics is Euler, which is set in
61 : : //! inciter::LuaParser::storeInputDeck()
62 : : template< class Physics, class Problem >
63 : : class MultiSpecies {
64 : :
65 : : private:
66 : : using eq = tag::multispecies;
67 : :
68 : : public:
69 : : //! Constructor
70 : 24 : explicit MultiSpecies() :
71 : : m_physics(),
72 : 24 : m_ncomp( g_inputdeck.get< tag::ncomp >() ),
73 : : m_nprim(nprim()),
74 : 24 : m_riemann( multispeciesRiemannSolver( g_inputdeck.get< tag::flux >() ) )
75 : : {
76 : : // associate boundary condition configurations with state functions
77 [ + - ][ + - ]: 312 : brigand::for_each< ctr::bclist::Keys >( ConfigBC( m_bc,
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + + ]
[ + + ][ - - ]
[ - - ]
78 : : // BC State functions
79 : : { dirichlet
80 : : , symmetry
81 : : , invalidBC // Outlet BC not implemented
82 : : , farfield
83 : : , extrapolate
84 : : , noslipwall },
85 : : // BC Gradient functions
86 : : { noOpGrad
87 : : , symmetryGrad
88 : : , noOpGrad
89 : : , noOpGrad
90 : : , noOpGrad
91 : : , noOpGrad }
92 : : ) );
93 : :
94 : : // EoS initialization
95 [ + - ]: 24 : initializeSpeciesEoS( m_mat_blk );
96 : 24 : }
97 : :
98 : : //! Find the number of primitive quantities required for this PDE system
99 : : //! \return The number of primitive quantities required to be stored for
100 : : //! this PDE system (zero for multi-species)
101 : 112 : std::size_t nprim() const { return 0; }
102 : :
103 : : //! Find the number of materials set up for this PDE system
104 : : //! \return The number of materials set up for this PDE system
105 : 0 : std::size_t nmat() const { return 1; }
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 : 88 : void numEquationDofs(std::vector< std::size_t >& numEqDof) const
111 : : {
112 : : // all equation-dofs initialized to ndofs
113 [ + + ]: 572 : for (std::size_t i=0; i<m_ncomp; ++i) {
114 : 484 : numEqDof.push_back(g_inputdeck.get< tag::ndof >());
115 : : }
116 : 88 : }
117 : :
118 : : //! Determine elements that lie inside the user-defined IC box
119 : : //! \param[in] geoElem Element geometry array
120 : : //! \param[in] nielem Number of internal elements
121 : : //! \param[in,out] inbox List of nodes at which box user ICs are set for
122 : : //! each IC box
123 : 88 : void IcBoxElems( const tk::Fields& geoElem,
124 : : std::size_t nielem,
125 : : std::vector< std::unordered_set< std::size_t > >& inbox ) const
126 : : {
127 : 88 : tk::BoxElems< eq >(geoElem, nielem, inbox);
128 : 88 : }
129 : :
130 : : //! Find how many 'stiff equations' in this PDE system
131 : : //! \return number of stiff equations. Zero for now, but will need to change
132 : : //! as chemical non-equilibrium is added
133 : 528 : std::size_t nstiffeq() const
134 : : {
135 : 528 : return 0;
136 : : }
137 : :
138 : : //! Find how many 'non-stiff equations' in this PDE system
139 : : //! \return number of non-stiff equations
140 : 176 : std::size_t nnonstiffeq() const
141 : : {
142 : 176 : return m_ncomp-nstiffeq();
143 : : }
144 : :
145 : : //! Locate the stiff equations.
146 : : //! \param[out] stiffEqIdx list with pointers to stiff equations. Empty
147 : : //! for now but will have to index to chemical non-equilibrium when added
148 : 0 : void setStiffEqIdx( std::vector< std::size_t >& stiffEqIdx ) const
149 : : {
150 : 0 : stiffEqIdx.resize(0);
151 : 0 : }
152 : :
153 : : //! Locate the nonstiff equations.
154 : : //! \param[out] nonStiffEqIdx list with pointers to nonstiff equations
155 : 0 : void setNonStiffEqIdx( std::vector< std::size_t >& nonStiffEqIdx ) const
156 : : {
157 : 0 : nonStiffEqIdx.resize(0);
158 : 0 : }
159 : :
160 : : //! Initialize the compressible flow equations, prepare for time integration
161 : : //! \param[in] L Block diagonal mass matrix
162 : : //! \param[in] inpoel Element-node connectivity
163 : : //! \param[in] coord Array of nodal coordinates
164 : : //! \param[in] inbox List of elements at which box user ICs are set for
165 : : //! each IC box
166 : : //! \param[in] elemblkid Element ids associated with mesh block ids where
167 : : //! user ICs are set
168 : : //! \param[in,out] unk Array of unknowns
169 : : //! \param[in] t Physical time
170 : : //! \param[in] nielem Number of internal elements
171 : 88 : void initialize( const tk::Fields& L,
172 : : const std::vector< std::size_t >& inpoel,
173 : : const tk::UnsMesh::Coords& coord,
174 : : const std::vector< std::unordered_set< std::size_t > >& inbox,
175 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
176 : : elemblkid,
177 : : tk::Fields& unk,
178 : : tk::real t,
179 : : const std::size_t nielem ) const
180 : : {
181 [ + - ][ + - ]: 88 : tk::initialize( m_ncomp, m_mat_blk, L, inpoel, coord,
182 : : Problem::initialize, unk, t, nielem );
183 : :
184 : 88 : const auto rdof = g_inputdeck.get< tag::rdof >();
185 : 88 : const auto& ic = g_inputdeck.get< tag::ic >();
186 : 88 : const auto& icbox = ic.get< tag::box >();
187 : 88 : const auto& icmbk = ic.get< tag::meshblock >();
188 : :
189 : : // Set initial conditions inside user-defined IC boxes and mesh blocks
190 [ + - ]: 176 : std::vector< tk::real > s(m_ncomp, 0.0);
191 [ + + ]: 9184 : for (std::size_t e=0; e<nielem; ++e) {
192 : : // inside user-defined box
193 [ + - ]: 9096 : if (!icbox.empty()) {
194 : 9096 : std::size_t bcnt = 0;
195 [ + + ]: 18192 : for (const auto& b : icbox) { // for all boxes
196 [ + - ][ + - ]: 9096 : if (inbox.size() > bcnt && inbox[bcnt].find(e) != inbox[bcnt].end())
[ + + ][ + + ]
197 : : {
198 [ + - ]: 9096 : std::vector< tk::real > box
199 : 4548 : { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
200 : 4548 : b.template get< tag::ymin >(), b.template get< tag::ymax >(),
201 : 4548 : b.template get< tag::zmin >(), b.template get< tag::zmax >() };
202 : 4548 : auto V_ex = (box[1]-box[0]) * (box[3]-box[2]) * (box[5]-box[4]);
203 [ + + ]: 29562 : for (std::size_t c=0; c<m_ncomp; ++c) {
204 : 25014 : auto mark = c*rdof;
205 [ + - ]: 25014 : s[c] = unk(e,mark);
206 : : // set high-order DOFs to zero
207 [ + + ]: 100056 : for (std::size_t i=1; i<rdof; ++i)
208 [ + - ]: 75042 : unk(e,mark+i) = 0.0;
209 : : }
210 [ + - ]: 4548 : initializeBox<ctr::boxList>( m_mat_blk, V_ex, t, b, s );
211 : : // store box-initialization in solution vector
212 [ + + ]: 29562 : for (std::size_t c=0; c<m_ncomp; ++c) {
213 : 25014 : auto mark = c*rdof;
214 [ + - ]: 25014 : unk(e,mark) = s[c];
215 : : }
216 : : }
217 : 9096 : ++bcnt;
218 : : }
219 : : }
220 : :
221 : : // inside user-specified mesh blocks
222 [ - + ]: 9096 : if (!icmbk.empty()) {
223 [ - - ]: 0 : for (const auto& b : icmbk) { // for all blocks
224 : 0 : auto blid = b.get< tag::blockid >();
225 : 0 : auto V_ex = b.get< tag::volume >();
226 [ - - ][ - - ]: 0 : if (elemblkid.find(blid) != elemblkid.end()) {
227 [ - - ]: 0 : const auto& elset = tk::cref_find(elemblkid, blid);
228 [ - - ][ - - ]: 0 : if (elset.find(e) != elset.end()) {
229 [ - - ]: 0 : initializeBox<ctr::meshblockList>( m_mat_blk, V_ex, t, b, s );
230 : : // store initialization in solution vector
231 [ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) {
232 : 0 : auto mark = c*rdof;
233 [ - - ]: 0 : unk(e,mark) = s[c];
234 : : }
235 : : }
236 : : }
237 : : }
238 : : }
239 : : }
240 : 88 : }
241 : :
242 : : //! Compute density constraint for a given material. No-op
243 : : // //! \param[in] nelem Number of elements
244 : : // //! \param[in] unk Array of unknowns
245 : : //! \param[out] densityConstr Density Constraint: rho/(rho0*det(g))
246 : 176 : void computeDensityConstr( std::size_t /*nelem*/,
247 : : tk::Fields& /*unk*/,
248 : : std::vector< tk::real >& densityConstr) const
249 : : {
250 : 176 : densityConstr.resize(0);
251 : 176 : }
252 : :
253 : : //! Compute the left hand side block-diagonal mass matrix
254 : : //! \param[in] geoElem Element geometry array
255 : : //! \param[in,out] l Block diagonal mass matrix
256 : 88 : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const {
257 : 88 : const auto ndof = g_inputdeck.get< tag::ndof >();
258 : 88 : tk::mass( m_ncomp, ndof, geoElem, l );
259 : 88 : }
260 : :
261 : : //! Update the interface cells to first order dofs. No-op.
262 : : // //! \param[in] unk Array of unknowns
263 : : // //! \param[in] nielem Number of internal elements
264 : : // //! \param[in,out] ndofel Array of dofs
265 : : // //! \param[in,out] interface Vector of interface marker
266 : 6600 : void updateInterfaceCells( tk::Fields& /*unk*/,
267 : : std::size_t /*nielem*/,
268 : : std::vector< std::size_t >& /*ndofel*/,
269 : 6600 : std::vector< std::size_t >& /*interface*/ ) const {}
270 : :
271 : : //! Update the primitives for this PDE system. No-op.
272 : : // //! \param[in] unk Array of unknowns
273 : : // //! \param[in] L The left hand side block-diagonal mass matrix
274 : : // //! \param[in] geoElem Element geometry array
275 : : // //! \param[in,out] prim Array of primitives
276 : : // //! \param[in] nielem Number of internal elements
277 : : // //! \param[in] ndofel Array of dofs
278 : 6688 : void updatePrimitives( const tk::Fields& /*unk*/,
279 : : const tk::Fields& /*L*/,
280 : : const tk::Fields& /*geoElem*/,
281 : : tk::Fields& /*prim*/,
282 : : std::size_t /*nielem*/,
283 : 6688 : std::vector< std::size_t >& /*ndofel*/ ) const {}
284 : :
285 : : //! Clean up the state of trace materials for this PDE system. No-op.
286 : : // //! \param[in] t Physical time
287 : : // //! \param[in] geoElem Element geometry array
288 : : // //! \param[in,out] unk Array of unknowns
289 : : // //! \param[in,out] prim Array of primitives
290 : : // //! \param[in] nielem Number of internal elements
291 : 6600 : void cleanTraceMaterial( tk::real /*t*/,
292 : : const tk::Fields& /*geoElem*/,
293 : : tk::Fields& /*unk*/,
294 : : tk::Fields& /*prim*/,
295 : 6600 : std::size_t /*nielem*/ ) const {}
296 : :
297 : : //! Reconstruct second-order solution from first-order
298 : : //! \param[in] geoElem Element geometry array
299 : : //! \param[in] fd Face connectivity and boundary conditions object
300 : : //! \param[in] esup Elements-surrounding-nodes connectivity
301 : : //! \param[in] inpoel Element-node connectivity
302 : : //! \param[in] coord Array of nodal coordinates
303 : : //! \param[in,out] U Solution vector at recent time step
304 : : // //! \param[in,out] P Vector of primitives at recent time step
305 : : //! \param[in] pref Indicator for p-adaptive algorithm
306 : : //! \param[in] ndofel Vector of local number of degrees of freedome
307 : 6600 : void reconstruct( tk::real,
308 : : const tk::Fields&,
309 : : const tk::Fields& geoElem,
310 : : const inciter::FaceData& fd,
311 : : const std::map< std::size_t, std::vector< std::size_t > >&
312 : : esup,
313 : : const std::vector< std::size_t >& inpoel,
314 : : const tk::UnsMesh::Coords& coord,
315 : : tk::Fields& U,
316 : : tk::Fields& /*P*/,
317 : : const bool pref,
318 : : const std::vector< std::size_t >& ndofel ) const
319 : : {
320 : 6600 : const auto rdof = g_inputdeck.get< tag::rdof >();
321 : 6600 : const auto ndof = g_inputdeck.get< tag::ndof >();
322 : :
323 : 6600 : bool is_p0p1(false);
324 [ + - ][ + - ]: 6600 : if (rdof == 4 && ndof == 1)
325 : 6600 : is_p0p1 = true;
326 : :
327 : 6600 : const auto nelem = fd.Esuel().size()/4;
328 : :
329 [ - + ][ - - ]: 6600 : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
330 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
331 : :
332 : : //----- reconstruction of conserved quantities -----
333 : : //--------------------------------------------------
334 : :
335 [ + + ]: 688800 : for (std::size_t e=0; e<nelem; ++e)
336 : : {
337 : 1364400 : std::vector< std::size_t > vars;
338 : : // check if element is marked as p0p1
339 [ - + ][ - - ]: 682200 : if ( (pref && ndofel[e] == 1) || is_p0p1 ) {
[ + - ][ + - ]
340 : : // 1. specify how many variables need to be reconstructed
341 [ + + ][ + - ]: 4434300 : for (std::size_t c=0; c<m_ncomp; ++c) vars.push_back(c);
342 : :
343 : : // 2. solve 3x3 least-squares system
344 : : // Reconstruct second-order dofs in Taylor space using nodal-stencils
345 [ + - ]: 682200 : tk::recoLeastSqExtStencil( rdof, e, esup, inpoel, geoElem, U, vars );
346 : :
347 : : // 3. transform reconstructed derivatives to Dubiner dofs
348 [ + - ]: 682200 : tk::transform_P0P1( rdof, e, inpoel, coord, U, vars );
349 : : }
350 : : }
351 : 6600 : }
352 : :
353 : : //! Limit second-order solution, and primitive quantities separately
354 : : // //! \param[in] pref Indicator for p-adaptive algorithm
355 : : //! \param[in] geoFace Face geometry array
356 : : //! \param[in] geoElem Element geometry array
357 : : //! \param[in] fd Face connectivity and boundary conditions object
358 : : //! \param[in] esup Elements-surrounding-nodes connectivity
359 : : //! \param[in] inpoel Element-node connectivity
360 : : //! \param[in] coord Array of nodal coordinates
361 : : //! \param[in] ndofel Vector of local number of degrees of freedome
362 : : // //! \param[in] gid Local->global node id map
363 : : // //! \param[in] bid Local chare-boundary node ids (value) associated to
364 : : // //! global node ids (key)
365 : : // //! \param[in] uNodalExtrm Chare-boundary nodal extrema for conservative
366 : : // //! variables
367 : : // //! \param[in] pNodalExtrm Chare-boundary nodal extrema for primitive
368 : : // //! variables
369 : : // //! \param[in] mtInv Inverse of Taylor mass matrix
370 : : //! \param[in,out] U Solution vector at recent time step
371 : : // //! \param[in,out] P Vector of primitives at recent time step
372 : : //! \param[in,out] shockmarker Vector of shock-marker values
373 : 6600 : void limit( [[maybe_unused]] tk::real,
374 : : const bool /*pref*/,
375 : : const tk::Fields& geoFace,
376 : : const tk::Fields& geoElem,
377 : : const inciter::FaceData& fd,
378 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
379 : : const std::vector< std::size_t >& inpoel,
380 : : const tk::UnsMesh::Coords& coord,
381 : : const std::vector< std::size_t >& ndofel,
382 : : const std::vector< std::size_t >& /*gid*/,
383 : : const std::unordered_map< std::size_t, std::size_t >& /*bid*/,
384 : : const std::vector< std::vector<tk::real> >& /*uNodalExtrm*/,
385 : : const std::vector< std::vector<tk::real> >& /*pNodalExtrm*/,
386 : : const std::vector< std::vector<tk::real> >& /*mtInv*/,
387 : : tk::Fields& U,
388 : : tk::Fields& /*P*/,
389 : : std::vector< std::size_t >& shockmarker ) const
390 : : {
391 : 6600 : const auto limiter = g_inputdeck.get< tag::limiter >();
392 : 6600 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
393 : 6600 : const auto rdof = g_inputdeck.get< tag::rdof >();
394 : 6600 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
395 : :
396 : : // limit vectors of conserved and primitive quantities
397 [ + - ][ + - ]: 6600 : if (limiter == ctr::LimiterType::VERTEXBASEDP1 && rdof == 4)
398 : : {
399 [ + - ]: 13200 : VertexBasedMultiSpecies_P1( esup, inpoel, ndofel, fd.Esuel().size()/4,
400 : 6600 : m_mat_blk, fd, geoFace, geoElem, coord, flux, solidx, U, nspec,
401 : : shockmarker );
402 : : }
403 [ - - ][ - - ]: 0 : else if (limiter == ctr::LimiterType::VERTEXBASEDP1 && rdof == 10)
404 : : {
405 [ - - ]: 0 : VertexBasedMultiSpecies_P2( esup, inpoel, ndofel, fd.Esuel().size()/4,
406 : 0 : m_mat_blk, fd, geoFace, geoElem, coord, flux, solidx, U, nspec,
407 : : shockmarker );
408 : : }
409 [ - - ]: 0 : else if (limiter != ctr::LimiterType::NOLIMITER)
410 : : {
411 [ - - ][ - - ]: 0 : Throw("Limiter type not configured for multispecies.");
[ - - ]
412 : : }
413 : 6600 : }
414 : :
415 : : //! Apply CPL to the conservative variable solution for this PDE system
416 : : // //! \param[in] prim Array of primitive variables
417 : : // //! \param[in] geoElem Element geometry array
418 : : // //! \param[in] inpoel Element-node connectivity
419 : : // //! \param[in] coord Array of nodal coordinates
420 : : // //! \param[in,out] unk Array of conservative variables
421 : : // //! \param[in] nielem Number of internal elements
422 : : //! \details This function applies CPL to obtain consistent dofs for
423 : : //! conservative quantities based on the limited primitive quantities.
424 : : //! No-op for now, but might need in the future, see appendix of paper.
425 : : //! See Pandare et al. (2023). On the Design of Stable,
426 : : //! Consistent, and Conservative High-Order Methods for Multi-Material
427 : : //! Hydrodynamics. J Comp Phys, 112313.
428 : 6600 : void CPL( const tk::Fields& /*prim*/,
429 : : const tk::Fields& /*geoElem*/,
430 : : const std::vector< std::size_t >& /*inpoel*/,
431 : : const tk::UnsMesh::Coords& /*coord*/,
432 : : tk::Fields& /*unk*/,
433 : 6600 : std::size_t /*nielem*/ ) const {}
434 : :
435 : : //! Return cell-average deformation gradient tensor. No-op.
436 : 0 : std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
437 : : const tk::Fields&,
438 : : std::size_t ) const
439 : 0 : { return {}; }
440 : :
441 : : //! Reset the high order solution for p-adaptive scheme
442 : : //! \param[in] fd Face connectivity and boundary conditions object
443 : : //! \param[in,out] unk Solution vector at recent time step
444 : : //! \param[in,out] prim Primitive vector at recent time step
445 : : //! \param[in] ndofel Vector of local number of degrees of freedome
446 : : //! \details This function reset the high order coefficient for p-adaptive
447 : : //! solution polynomials. Unlike compflow class, the high order of fv
448 : : //! solution will not be reset since p0p1 is the base scheme for
449 : : //! multi-species p-adaptive DG method.
450 : 0 : void resetAdapSol( const inciter::FaceData& fd,
451 : : tk::Fields& unk,
452 : : tk::Fields& prim,
453 : : const std::vector< std::size_t >& ndofel ) const
454 : : {
455 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
456 : 0 : const auto ncomp = unk.nprop() / rdof;
457 : 0 : const auto nprim = prim.nprop() / rdof;
458 : :
459 [ - - ]: 0 : for(std::size_t e = 0; e < fd.Esuel().size()/4; e++)
460 : : {
461 [ - - ]: 0 : if(ndofel[e] < 10)
462 : : {
463 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c)
464 : : {
465 : 0 : auto mark = c*rdof;
466 : 0 : unk(e, mark+4) = 0.0;
467 : 0 : unk(e, mark+5) = 0.0;
468 : 0 : unk(e, mark+6) = 0.0;
469 : 0 : unk(e, mark+7) = 0.0;
470 : 0 : unk(e, mark+8) = 0.0;
471 : 0 : unk(e, mark+9) = 0.0;
472 : : }
473 [ - - ]: 0 : for (std::size_t c=0; c<nprim; ++c)
474 : : {
475 : 0 : auto mark = c*rdof;
476 : 0 : prim(e, mark+4) = 0.0;
477 : 0 : prim(e, mark+5) = 0.0;
478 : 0 : prim(e, mark+6) = 0.0;
479 : 0 : prim(e, mark+7) = 0.0;
480 : 0 : prim(e, mark+8) = 0.0;
481 : 0 : prim(e, mark+9) = 0.0;
482 : : }
483 : : }
484 : : }
485 : 0 : }
486 : :
487 : : //! Compute right hand side
488 : : //! \param[in] t Physical time
489 : : //! \param[in] pref Indicator for p-adaptive algorithm
490 : : //! \param[in] geoFace Face geometry array
491 : : //! \param[in] geoElem Element geometry array
492 : : //! \param[in] fd Face connectivity and boundary conditions object
493 : : //! \param[in] inpoel Element-node connectivity
494 : : //! \param[in] coord Array of nodal coordinates
495 : : //! \param[in] U Solution vector at recent time step
496 : : //! \param[in] P Primitive vector at recent time step
497 : : //! \param[in] ndofel Vector of local number of degrees of freedom
498 : : //! \param[in] dt Delta time
499 : : //! \param[in,out] R Right-hand side vector computed
500 : 6600 : void rhs( tk::real t,
501 : : const bool pref,
502 : : const tk::Fields& geoFace,
503 : : const tk::Fields& geoElem,
504 : : const inciter::FaceData& fd,
505 : : const std::vector< std::size_t >& inpoel,
506 : : const std::vector< std::unordered_set< std::size_t > >&,
507 : : const tk::UnsMesh::Coords& coord,
508 : : const tk::Fields& U,
509 : : const tk::Fields& P,
510 : : const std::vector< std::size_t >& ndofel,
511 : : const tk::real dt,
512 : : tk::Fields& R ) const
513 : : {
514 : 6600 : const auto ndof = g_inputdeck.get< tag::ndof >();
515 : 6600 : const auto rdof = g_inputdeck.get< tag::rdof >();
516 : 6600 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
517 : :
518 : 6600 : const auto nelem = fd.Esuel().size()/4;
519 : :
520 [ - + ][ - - ]: 6600 : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
521 : : "vector and primitive vector at recent time step incorrect" );
522 [ - + ][ - - ]: 6600 : Assert( U.nunk() == R.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
523 : : "vector and right-hand side at recent time step incorrect" );
524 [ - + ][ - - ]: 6600 : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
525 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
526 [ - + ][ - - ]: 6600 : Assert( P.nprop() == rdof*m_nprim, "Number of components in primitive "
[ - - ][ - - ]
[ - - ]
527 : : "vector must equal "+ std::to_string(rdof*m_nprim) );
528 [ - + ][ - - ]: 6600 : Assert( R.nprop() == ndof*m_ncomp, "Number of components in right-hand "
[ - - ][ - - ]
[ - - ]
529 : : "side vector must equal "+ std::to_string(ndof*m_ncomp) );
530 [ - + ][ - - ]: 6600 : Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
[ - - ][ - - ]
531 : : "Mismatch in inpofa size" );
532 : :
533 : : // set rhs to zero
534 [ + - ]: 6600 : R.fill(0.0);
535 : :
536 : : // empty vector for non-conservative terms. This vector is unused for
537 : : // multi-species flow since, there are no non-conservative terms
538 : : // in the system of PDEs.
539 : 13200 : std::vector< std::vector< tk::real > > riemannDeriv;
540 : :
541 : 13200 : std::vector< std::vector< tk::real > > vriem;
542 : 13200 : std::vector< std::vector< tk::real > > riemannLoc;
543 : :
544 : : // configure a no-op lambda for prescribed velocity
545 : 1752300 : auto velfn = []( ncomp_t, tk::real, tk::real, tk::real, tk::real ){
546 : 1752300 : return tk::VelFn::result_type(); };
547 : :
548 : : // compute internal surface flux integrals
549 [ + - ]: 13200 : tk::surfInt( pref, 1, m_mat_blk, t, ndof, rdof, inpoel, solidx,
550 [ + - ]: 6600 : coord, fd, geoFace, geoElem, m_riemann, velfn, U, P, ndofel,
551 : : dt, R, riemannDeriv );
552 : :
553 : : // compute optional source term
554 [ + - ][ + - ]: 6600 : tk::srcInt( m_mat_blk, t, ndof, fd.Esuel().size()/4, inpoel,
555 : : coord, geoElem, Problem::src, ndofel, R );
556 : :
557 [ - + ]: 6600 : if(ndof > 1)
558 : : // compute volume integrals
559 [ - - ][ - - ]: 0 : tk::volInt( 1, t, m_mat_blk, ndof, rdof, nelem, inpoel, coord, geoElem,
[ - - ]
560 : : flux, velfn, U, P, ndofel, R );
561 : :
562 : : // compute boundary surface flux integrals
563 [ + + ]: 46200 : for (const auto& b : m_bc)
564 [ + - ][ + - ]: 79200 : tk::bndSurfInt( pref, 1, m_mat_blk, ndof, rdof, std::get<0>(b), fd,
565 : 39600 : geoFace, geoElem, inpoel, coord, t, m_riemann, velfn,
566 : 39600 : std::get<1>(b), U, P, ndofel, R, riemannDeriv );
567 : :
568 : : // compute external (energy) sources
569 : : //m_physics.physSrc(nspec, t, geoElem, {}, R, {});
570 : 6600 : }
571 : :
572 : : //! Evaluate the adaptive indicator and mark the ndof for each element
573 : : //! \param[in] nunk Number of unknowns
574 : : //! \param[in] coord Array of nodal coordinates
575 : : //! \param[in] inpoel Element-node connectivity
576 : : //! \param[in] fd Face connectivity and boundary conditions object
577 : : //! \param[in] unk Array of unknowns
578 : : // //! \param[in] prim Array of primitive quantities
579 : : //! \param[in] indicator p-refinement indicator type
580 : : //! \param[in] ndof Number of degrees of freedom in the solution
581 : : //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
582 : : //! \param[in] tolref Tolerance for p-refinement
583 : : //! \param[in,out] ndofel Vector of local number of degrees of freedome
584 : 0 : void eval_ndof( std::size_t nunk,
585 : : [[maybe_unused]] const tk::UnsMesh::Coords& coord,
586 : : [[maybe_unused]] const std::vector< std::size_t >& inpoel,
587 : : const inciter::FaceData& fd,
588 : : const tk::Fields& unk,
589 : : const tk::Fields& /*prim*/,
590 : : inciter::ctr::PrefIndicatorType indicator,
591 : : std::size_t ndof,
592 : : std::size_t ndofmax,
593 : : tk::real tolref,
594 : : std::vector< std::size_t >& ndofel ) const
595 : : {
596 : 0 : const auto& esuel = fd.Esuel();
597 : :
598 [ - - ]: 0 : if(indicator == inciter::ctr::PrefIndicatorType::SPECTRAL_DECAY)
599 : 0 : spectral_decay(1, nunk, esuel, unk, ndof, ndofmax, tolref, ndofel);
600 : : else
601 [ - - ][ - - ]: 0 : Throw( "No such adaptive indicator type" );
[ - - ]
602 : 0 : }
603 : :
604 : : //! Compute the minimum time step size
605 : : //! \param[in] fd Face connectivity and boundary conditions object
606 : : //! \param[in] geoFace Face geometry array
607 : : //! \param[in] geoElem Element geometry array
608 : : // //! \param[in] ndofel Vector of local number of degrees of freedom
609 : : //! \param[in] U Solution vector at recent time step
610 : : //! \param[in] P Vector of primitive quantities at recent time step
611 : : //! \param[in] nielem Number of internal elements
612 : : //! \return Minimum time step size
613 : : //! \details The allowable dt is calculated by looking at the maximum
614 : : //! wave-speed in elements surrounding each face, times the area of that
615 : : //! face. Once the maximum of this quantity over the mesh is determined,
616 : : //! the volume of each cell is divided by this quantity. A minimum of this
617 : : //! ratio is found over the entire mesh, which gives the allowable dt.
618 : 0 : tk::real dt( const std::array< std::vector< tk::real >, 3 >&,
619 : : const std::vector< std::size_t >&,
620 : : const inciter::FaceData& fd,
621 : : const tk::Fields& geoFace,
622 : : const tk::Fields& geoElem,
623 : : const std::vector< std::size_t >& /*ndofel*/,
624 : : const tk::Fields& U,
625 : : const tk::Fields& P,
626 : : const std::size_t nielem ) const
627 : : {
628 : 0 : const auto ndof = g_inputdeck.get< tag::ndof >();
629 : 0 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
630 : :
631 : 0 : auto mindt = timeStepSizeMultiSpecies( m_mat_blk, fd.Esuf(), geoFace,
632 : : geoElem, nielem, nspec, U, P);
633 : :
634 : : //if (viscous)
635 : : // mindt = std::min(mindt, timeStepSizeViscousFV(geoElem, nielem, nspec, U));
636 : : //mindt = std::min(mindt, m_physics.dtRestriction(geoElem, nielem, {}));
637 : :
638 : 0 : tk::real dgp = 0.0;
639 [ - - ]: 0 : if (ndof == 4)
640 : : {
641 : 0 : dgp = 1.0;
642 : : }
643 [ - - ]: 0 : else if (ndof == 10)
644 : : {
645 : 0 : dgp = 2.0;
646 : : }
647 : :
648 : : // Scale smallest dt with CFL coefficient and the CFL is scaled by (2*p+1)
649 : : // where p is the order of the DG polynomial by linear stability theory.
650 : 0 : mindt /= (2.0*dgp + 1.0);
651 : 0 : return mindt;
652 : : }
653 : :
654 : : //! Compute stiff terms for a single element. No-op until chem sources added
655 : : // //! \param[in] e Element number
656 : : // //! \param[in] geoElem Element geometry array
657 : : // //! \param[in] inpoel Element-node connectivity
658 : : // //! \param[in] coord Array of nodal coordinates
659 : : // //! \param[in] U Solution vector at recent time step
660 : : // //! \param[in] P Primitive vector at recent time step
661 : : // //! \param[in] ndofel Vector of local number of degrees of freedom
662 : : // //! \param[in,out] R Right-hand side vector computed
663 : 0 : void stiff_rhs( std::size_t /*e*/,
664 : : const tk::Fields& /*geoElem*/,
665 : : const std::vector< std::size_t >& /*inpoel*/,
666 : : const tk::UnsMesh::Coords& /*coord*/,
667 : : const tk::Fields& /*U*/,
668 : : const tk::Fields& /*P*/,
669 : : const std::vector< std::size_t >& /*ndofel*/,
670 : 0 : tk::Fields& /*R*/ ) const {}
671 : :
672 : : //! Extract the velocity field at cell nodes. Currently unused.
673 : : // //! \param[in] U Solution vector at recent time step
674 : : // //! \param[in] N Element node indices
675 : : //! \return Array of the four values of the velocity field
676 : : std::array< std::array< tk::real, 4 >, 3 >
677 : : velocity( const tk::Fields& /*U*/,
678 : : const std::array< std::vector< tk::real >, 3 >&,
679 : : const std::array< std::size_t, 4 >& /*N*/ ) const
680 : : {
681 : : std::array< std::array< tk::real, 4 >, 3 > v;
682 : : return v;
683 : : }
684 : :
685 : : //! Return a map that associates user-specified strings to functions
686 : : //! \return Map that associates user-specified strings to functions that
687 : : //! compute relevant quantities to be output to file
688 : 352 : std::map< std::string, tk::GetVarFn > OutVarFn() const
689 : 352 : { return MultiSpeciesOutVarFn(); }
690 : :
691 : : //! Return analytic field names to be output to file
692 : : //! \return Vector of strings labelling analytic fields output in file
693 : 0 : std::vector< std::string > analyticFieldNames() const {
694 : 0 : auto nspec = g_inputdeck.get< eq, tag::nspec >();
695 : :
696 : 0 : return MultiSpeciesFieldNames(nspec);
697 : : }
698 : :
699 : : //! Return time history field names to be output to file
700 : : //! \return Vector of strings labelling time history fields output in file
701 : 0 : std::vector< std::string > histNames() const {
702 : 0 : return MultiSpeciesHistNames();
703 : : }
704 : :
705 : : //! Return surface field output going to file
706 : : std::vector< std::vector< tk::real > >
707 : 0 : surfOutput( const std::map< int, std::vector< std::size_t > >&,
708 : : tk::Fields& ) const
709 : : {
710 : 0 : std::vector< std::vector< tk::real > > s; // punt for now
711 : 0 : return s;
712 : : }
713 : :
714 : : //! Return time history field output evaluated at time history points
715 : : //! \param[in] h History point data
716 : : //! \param[in] inpoel Element-node connectivity
717 : : //! \param[in] coord Array of nodal coordinates
718 : : //! \param[in] U Array of unknowns
719 : : // //! \param[in] P Array of primitive quantities
720 : : //! \return Vector of time history output of bulk flow quantities (density,
721 : : //! velocity, total energy, and pressure) evaluated at time history points
722 : : std::vector< std::vector< tk::real > >
723 : 0 : histOutput( const std::vector< HistData >& h,
724 : : const std::vector< std::size_t >& inpoel,
725 : : const tk::UnsMesh::Coords& coord,
726 : : const tk::Fields& U,
727 : : const tk::Fields& /*P*/ ) const
728 : : {
729 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
730 : 0 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
731 : :
732 : 0 : const auto& x = coord[0];
733 : 0 : const auto& y = coord[1];
734 : 0 : const auto& z = coord[2];
735 : :
736 [ - - ]: 0 : std::vector< std::vector< tk::real > > Up(h.size());
737 : :
738 : 0 : std::size_t j = 0;
739 [ - - ]: 0 : for (const auto& p : h) {
740 : 0 : auto e = p.get< tag::elem >();
741 : 0 : auto chp = p.get< tag::coord >();
742 : :
743 : : // Evaluate inverse Jacobian
744 : 0 : std::array< std::array< tk::real, 3>, 4 > cp{{
745 : : {{ x[inpoel[4*e ]], y[inpoel[4*e ]], z[inpoel[4*e ]] }},
746 : 0 : {{ x[inpoel[4*e+1]], y[inpoel[4*e+1]], z[inpoel[4*e+1]] }},
747 : 0 : {{ x[inpoel[4*e+2]], y[inpoel[4*e+2]], z[inpoel[4*e+2]] }},
748 : 0 : {{ x[inpoel[4*e+3]], y[inpoel[4*e+3]], z[inpoel[4*e+3]] }} }};
749 : 0 : auto J = tk::inverseJacobian( cp[0], cp[1], cp[2], cp[3] );
750 : :
751 : : // evaluate solution at history-point
752 : 0 : std::array< tk::real, 3 > dc{{chp[0]-cp[0][0], chp[1]-cp[0][1],
753 : 0 : chp[2]-cp[0][2]}};
754 [ - - ]: 0 : auto B = tk::eval_basis(rdof, tk::dot(J[0],dc), tk::dot(J[1],dc),
755 : 0 : tk::dot(J[2],dc));
756 [ - - ]: 0 : auto uhp = eval_state(m_ncomp, rdof, rdof, e, U, B);
757 : :
758 : : // store solution in history output vector
759 [ - - ]: 0 : Up[j].resize(6+nspec, 0.0);
760 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k) {
761 : 0 : Up[j][0] += uhp[multispecies::densityIdx(nspec,k)];
762 : : }
763 : 0 : Up[j][1] = uhp[multispecies::momentumIdx(nspec,0)]/Up[j][0];
764 : 0 : Up[j][2] = uhp[multispecies::momentumIdx(nspec,1)]/Up[j][0];
765 : 0 : Up[j][3] = uhp[multispecies::momentumIdx(nspec,2)]/Up[j][0];
766 : 0 : Up[j][4] = uhp[multispecies::energyIdx(nspec,0)];
767 [ - - ]: 0 : Up[j][5] = m_mat_blk[0].compute< EOS::pressure >( Up[j][0], Up[j][1],
768 : : Up[j][2], Up[j][3], Up[j][4]);
769 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k) {
770 : 0 : Up[j][6+k] = uhp[multispecies::densityIdx(nspec,k)]/Up[j][0];
771 : : }
772 : 0 : ++j;
773 : : }
774 : :
775 : 0 : return Up;
776 : : }
777 : :
778 : : //! Return names of integral variables to be output to diagnostics file
779 : : //! \return Vector of strings labelling integral variables output
780 : 6 : std::vector< std::string > names() const
781 : : {
782 : 6 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
783 : 6 : return MultiSpeciesDiagNames(nspec);
784 : : }
785 : :
786 : : //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
787 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
788 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
789 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
790 : : //! \param[in] t Physical time at which to evaluate the analytic solution
791 : : //! \return Vector of analytic solution at given location and time
792 : : std::vector< tk::real >
793 : 0 : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
794 : 0 : { return Problem::analyticSolution( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
795 : :
796 : : //! Return analytic solution for conserved variables
797 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
798 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
799 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
800 : : //! \param[in] t Physical time at which to evaluate the analytic solution
801 : : //! \return Vector of analytic solution at given location and time
802 : : std::vector< tk::real >
803 : 227400 : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
804 : 227400 : { return Problem::initialize( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
805 : :
806 : : //! Return cell-averaged specific total energy for an element
807 : : //! \param[in] e Element id for which total energy is required
808 : : //! \param[in] unk Vector of conserved quantities
809 : : //! \return Cell-averaged specific total energy for given element
810 : 227400 : tk::real sp_totalenergy(std::size_t e, const tk::Fields& unk) const
811 : : {
812 : 227400 : const auto rdof = g_inputdeck.get< tag::rdof >();
813 : 227400 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
814 : :
815 : 227400 : return unk(e, multispecies::energyDofIdx(nspec,0,rdof,0));
816 : : }
817 : :
818 : : private:
819 : : //! Physics policy
820 : : const Physics m_physics;
821 : : //! Number of components in this PDE system
822 : : const ncomp_t m_ncomp;
823 : : //! Number of primitive quantities stored in this PDE system
824 : : const ncomp_t m_nprim;
825 : : //! Riemann solver
826 : : tk::RiemannFluxFn m_riemann;
827 : : //! BC configuration
828 : : BCStateFn m_bc;
829 : : //! EOS material block
830 : : std::vector< EOS > m_mat_blk;
831 : :
832 : : //! Evaluate conservative part of physical flux function for this PDE system
833 : : //! \param[in] ncomp Number of scalar components in this PDE system
834 : : //! \param[in] ugp Numerical solution at the Gauss point at which to
835 : : //! evaluate the flux
836 : : //! \return Flux vectors for all components in this PDE system
837 : : //! \note The function signature must follow tk::FluxFn
838 : : static tk::FluxFn::result_type
839 : 0 : flux( [[maybe_unused]] ncomp_t ncomp,
840 : : const std::vector< EOS >& mat_blk,
841 : : const std::vector< tk::real >& ugp,
842 : : const std::vector< std::array< tk::real, 3 > >& )
843 : : {
844 [ - - ][ - - ]: 0 : Assert( ugp.size() == ncomp, "Size mismatch" );
[ - - ][ - - ]
845 : :
846 : 0 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
847 : :
848 [ - - ]: 0 : std::vector< std::array< tk::real, 3 > > fl( ugp.size() );
849 : :
850 : 0 : tk::real rhob(0.0);
851 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k)
852 : 0 : rhob += ugp[multispecies::densityIdx(nspec, k)];
853 : :
854 : 0 : std::array< tk::real, 3 > u{{
855 : 0 : ugp[multispecies::momentumIdx(nspec,0)] / rhob,
856 : 0 : ugp[multispecies::momentumIdx(nspec,1)] / rhob,
857 : 0 : ugp[multispecies::momentumIdx(nspec,2)] / rhob }};
858 : 0 : auto rhoE0 = ugp[multispecies::energyIdx(nspec,0)];
859 [ - - ]: 0 : auto p = mat_blk[0].compute< EOS::pressure >( rhob, u[0], u[1], u[2],
860 : : rhoE0 );
861 : :
862 : : // density flux
863 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k) {
864 : 0 : auto idx = multispecies::densityIdx(nspec, k);
865 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
866 : 0 : fl[idx][j] = ugp[idx] * u[j];
867 : : }
868 : : }
869 : :
870 : : // momentum flux
871 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
872 : 0 : auto idx = multispecies::momentumIdx(nspec,i);
873 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
874 : 0 : fl[idx][j] = ugp[idx] * u[j];
875 [ - - ]: 0 : if (i == j) fl[idx][j] += p;
876 : : }
877 : : }
878 : :
879 : : // energy flux
880 : 0 : auto idx = multispecies::energyIdx(nspec,0);
881 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
882 : 0 : fl[idx][j] = u[j] * (ugp[idx] + p);
883 : : }
884 : :
885 : 0 : return fl;
886 : : }
887 : :
888 : : //! \brief Boundary state function providing the left and right state of a
889 : : //! face at Dirichlet boundaries
890 : : //! \param[in] ncomp Number of scalar components in this PDE system
891 : : //! \param[in] mat_blk EOS material block
892 : : //! \param[in] ul Left (domain-internal) state
893 : : //! \param[in] x X-coordinate at which to compute the states
894 : : //! \param[in] y Y-coordinate at which to compute the states
895 : : //! \param[in] z Z-coordinate at which to compute the states
896 : : //! \param[in] t Physical time
897 : : //! \return Left and right states for all scalar components in this PDE
898 : : //! system
899 : : //! \note The function signature must follow tk::StateFn.
900 : : static tk::StateFn::result_type
901 : 2250 : dirichlet( ncomp_t ncomp,
902 : : const std::vector< EOS >& mat_blk,
903 : : const std::vector< tk::real >& ul, tk::real x, tk::real y,
904 : : tk::real z, tk::real t, const std::array< tk::real, 3 >& )
905 : : {
906 : 2250 : return {{ ul, Problem::initialize( ncomp, mat_blk, x, y, z, t ) }};
907 : : }
908 : :
909 : : // Other boundary condition types that do not depend on "Problem" should be
910 : : // added in BCFunctions.hpp
911 : : };
912 : :
913 : : } // dg::
914 : :
915 : : } // inciter::
916 : :
917 : : #endif // DGMultiSpecies_h
|
