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/Surface.hpp"
33 : : #include "Integrate/Boundary.hpp"
34 : : #include "Integrate/Volume.hpp"
35 : : #include "Integrate/Source.hpp"
36 : : #include "Integrate/SolidTerms.hpp"
37 : : #include "RiemannChoice.hpp"
38 : : #include "MultiSpecies/MultiSpeciesIndexing.hpp"
39 : : #include "Reconstruction.hpp"
40 : : #include "Limiter.hpp"
41 : : #include "Problem/FieldOutput.hpp"
42 : : #include "Problem/BoxInitialization.hpp"
43 : : #include "PrefIndicator.hpp"
44 : : #include "MultiSpecies/BCFunctions.hpp"
45 : : #include "MultiSpecies/MiscMultiSpeciesFns.hpp"
46 : : #include "EoS/GetMatProp.hpp"
47 : : #include "Mixture/Mixture.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 [ + - ][ + - ]: 360 : 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 : : , symmetry }, // Slip equivalent to symmetry without mesh motion
86 : : // BC Gradient functions
87 : : { noOpGrad
88 : : , symmetryGrad
89 : : , noOpGrad
90 : : , noOpGrad
91 : : , noOpGrad
92 : : , noOpGrad
93 : : , symmetryGrad }
94 : : ) );
95 : :
96 : : // EoS initialization
97 [ + - ]: 24 : initializeSpeciesEoS( m_mat_blk );
98 : 24 : }
99 : :
100 : : //! Find the number of primitive quantities required for this PDE system
101 : : //! \return The number of primitive quantities required to be stored for
102 : : //! this PDE system
103 : 112 : std::size_t nprim() const
104 : : {
105 : : // mixture temperature
106 : 112 : std::size_t np(1);
107 : :
108 : 112 : return np;
109 : : }
110 : :
111 : : //! Find the number of materials set up for this PDE system
112 : : //! \return The number of materials set up for this PDE system
113 : 0 : std::size_t nmat() const { return 1; }
114 : :
115 : : //! Assign number of DOFs per equation in the PDE system
116 : : //! \param[in,out] numEqDof Array storing number of Dofs for each PDE
117 : : //! equation
118 : 88 : void numEquationDofs(std::vector< std::size_t >& numEqDof) const
119 : : {
120 : : // all equation-dofs initialized to ndofs
121 [ + + ]: 572 : for (std::size_t i=0; i<m_ncomp; ++i) {
122 : 484 : numEqDof.push_back(g_inputdeck.get< tag::ndof >());
123 : : }
124 : 88 : }
125 : :
126 : : //! Determine elements that lie inside the user-defined IC box
127 : : //! \param[in] geoElem Element geometry array
128 : : //! \param[in] nielem Number of internal elements
129 : : //! \param[in,out] inbox List of nodes at which box user ICs are set for
130 : : //! each IC box
131 : 88 : void IcBoxElems( const tk::Fields& geoElem,
132 : : std::size_t nielem,
133 : : std::vector< std::unordered_set< std::size_t > >& inbox ) const
134 : : {
135 : 88 : tk::BoxElems< eq >(geoElem, nielem, inbox);
136 : 88 : }
137 : :
138 : : //! Find how many 'stiff equations' in this PDE system
139 : : //! \return number of stiff equations. Zero for now, but will need to change
140 : : //! as chemical non-equilibrium is added
141 : 528 : std::size_t nstiffeq() const
142 : : {
143 : 528 : return 0;
144 : : }
145 : :
146 : : //! Find how many 'non-stiff equations' in this PDE system
147 : : //! \return number of non-stiff equations
148 : 176 : std::size_t nnonstiffeq() const
149 : : {
150 : 176 : return m_ncomp-nstiffeq();
151 : : }
152 : :
153 : : //! Locate the stiff equations.
154 : : //! \param[out] stiffEqIdx list with pointers to stiff equations. Empty
155 : : //! for now but will have to index to chemical non-equilibrium when added
156 : 0 : void setStiffEqIdx( std::vector< std::size_t >& stiffEqIdx ) const
157 : : {
158 : 0 : stiffEqIdx.resize(0);
159 : 0 : }
160 : :
161 : : //! Locate the nonstiff equations.
162 : : //! \param[out] nonStiffEqIdx list with pointers to nonstiff equations
163 : 0 : void setNonStiffEqIdx( std::vector< std::size_t >& nonStiffEqIdx ) const
164 : : {
165 : 0 : nonStiffEqIdx.resize(0);
166 : 0 : }
167 : :
168 : : //! Initialize the compressible flow equations, prepare for time integration
169 : : //! \param[in] geoElem Element geometry array
170 : : //! \param[in] inpoel Element-node connectivity
171 : : //! \param[in] coord Array of nodal coordinates
172 : : //! \param[in] inbox List of elements at which box user ICs are set for
173 : : //! each IC box
174 : : //! \param[in] elemblkid Element ids associated with mesh block ids where
175 : : //! user ICs are set
176 : : //! \param[in,out] unk Array of unknowns
177 : : //! \param[in] t Physical time
178 : : //! \param[in] nielem Number of internal elements
179 : 88 : void initialize( const tk::Fields& geoElem,
180 : : const std::vector< std::size_t >& inpoel,
181 : : const tk::UnsMesh::Coords& coord,
182 : : const std::vector< std::unordered_set< std::size_t > >& inbox,
183 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
184 : : elemblkid,
185 : : tk::Fields& unk,
186 : : tk::real t,
187 : : const std::size_t nielem ) const
188 : : {
189 [ + - ][ + - ]: 88 : tk::initialize( m_ncomp, m_mat_blk, geoElem, inpoel, coord,
190 : : Problem::initialize, unk, t, nielem );
191 : :
192 : 88 : const auto rdof = g_inputdeck.get< tag::rdof >();
193 : 88 : const auto& ic = g_inputdeck.get< tag::ic >();
194 : 88 : const auto& icbox = ic.get< tag::box >();
195 : 88 : const auto& icmbk = ic.get< tag::meshblock >();
196 : :
197 : : // Set initial conditions inside user-defined IC boxes and mesh blocks
198 [ + - ]: 176 : std::vector< tk::real > s(m_ncomp, 0.0);
199 [ + + ]: 9184 : for (std::size_t e=0; e<nielem; ++e) {
200 : : // inside user-defined box
201 [ + - ]: 9096 : if (!icbox.empty()) {
202 : 9096 : std::size_t bcnt = 0;
203 [ + + ]: 18192 : for (const auto& b : icbox) { // for all boxes
204 [ + - ][ + - ]: 9096 : if (inbox.size() > bcnt && inbox[bcnt].find(e) != inbox[bcnt].end())
[ + + ][ + + ]
205 : : {
206 [ + - ]: 9096 : std::vector< tk::real > box
207 : 4548 : { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
208 : 4548 : b.template get< tag::ymin >(), b.template get< tag::ymax >(),
209 : 4548 : b.template get< tag::zmin >(), b.template get< tag::zmax >() };
210 : 4548 : auto V_ex = (box[1]-box[0]) * (box[3]-box[2]) * (box[5]-box[4]);
211 [ + + ]: 29562 : for (std::size_t c=0; c<m_ncomp; ++c) {
212 : 25014 : auto mark = c*rdof;
213 [ + - ]: 25014 : s[c] = unk(e,mark);
214 : : // set high-order DOFs to zero
215 [ + + ]: 100056 : for (std::size_t i=1; i<rdof; ++i)
216 [ + - ]: 75042 : unk(e,mark+i) = 0.0;
217 : : }
218 [ + - ]: 4548 : initializeBox<ctr::boxList>( m_mat_blk, V_ex, t, b, s );
219 : : // store box-initialization in solution vector
220 [ + + ]: 29562 : for (std::size_t c=0; c<m_ncomp; ++c) {
221 : 25014 : auto mark = c*rdof;
222 [ + - ]: 25014 : unk(e,mark) = s[c];
223 : : }
224 : : }
225 : 9096 : ++bcnt;
226 : : }
227 : : }
228 : :
229 : : // inside user-specified mesh blocks
230 [ - + ]: 9096 : if (!icmbk.empty()) {
231 [ - - ]: 0 : for (const auto& b : icmbk) { // for all blocks
232 : 0 : auto blid = b.get< tag::blockid >();
233 : 0 : auto V_ex = b.get< tag::volume >();
234 [ - - ][ - - ]: 0 : if (elemblkid.find(blid) != elemblkid.end()) {
235 [ - - ]: 0 : const auto& elset = tk::cref_find(elemblkid, blid);
236 [ - - ][ - - ]: 0 : if (elset.find(e) != elset.end()) {
237 [ - - ]: 0 : initializeBox<ctr::meshblockList>( m_mat_blk, V_ex, t, b, s );
238 : : // store initialization in solution vector
239 [ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) {
240 : 0 : auto mark = c*rdof;
241 [ - - ]: 0 : unk(e,mark) = s[c];
242 : : }
243 : : }
244 : : }
245 : : }
246 : : }
247 : : }
248 : 88 : }
249 : :
250 : : //! Compute density constraint for a given material. No-op
251 : : // //! \param[in] nelem Number of elements
252 : : // //! \param[in] unk Array of unknowns
253 : : //! \param[out] densityConstr Density Constraint: rho/(rho0*det(g))
254 : 176 : void computeDensityConstr( std::size_t /*nelem*/,
255 : : tk::Fields& /*unk*/,
256 : : std::vector< tk::real >& densityConstr) const
257 : : {
258 : 176 : densityConstr.resize(0);
259 : 176 : }
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.
272 : : //! \param[in] unk Array of unknowns
273 : : //! \param[in] geoElem Element geometry array
274 : : //! \param[in,out] prim Array of primitives
275 : : //! \param[in] nielem Number of internal elements
276 : : //! \param[in] ndofel Array of dofs
277 : : //! \details This function computes and stores the dofs for primitive
278 : : //! quantities, which are required for obtaining reconstructed states used
279 : : //! in the Riemann solver. See for eg. /PDE/Riemann/AUSMMultiSpecies.hpp,
280 : : //! where temperature is independently reconstructed.
281 : 6688 : void updatePrimitives( const tk::Fields& unk,
282 : : const tk::Fields& geoElem,
283 : : tk::Fields& prim,
284 : : std::size_t nielem,
285 : : const std::vector< std::size_t >& ndofel ) const
286 : : {
287 : 6688 : const auto rdof = g_inputdeck.get< tag::rdof >();
288 : 6688 : const auto ndof = g_inputdeck.get< tag::ndof >();
289 : 6688 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
290 : :
291 [ - + ][ - - ]: 6688 : Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
292 : : "vector and primitive vector at recent time step incorrect" );
293 [ - + ][ - - ]: 6688 : Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
294 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
295 [ - + ][ - - ]: 6688 : Assert( prim.nprop() == rdof*m_nprim, "Number of components in vector of "
[ - - ][ - - ]
[ - - ]
296 : : "primitive quantities must equal "+ std::to_string(rdof*m_nprim) );
297 : :
298 : 6688 : auto mass_m = tk::massMatrixDubiner();
299 : :
300 [ + + ]: 697984 : for (std::size_t e=0; e<nielem; ++e)
301 : : {
302 [ + - ]: 1382592 : std::vector< tk::real > R(m_nprim*ndof, 0.0);
303 : :
304 [ + - ]: 691296 : auto ng = tk::NGvol(ndof);
305 : :
306 : : // arrays for quadrature points
307 : 1382592 : std::array< std::vector< tk::real >, 3 > coordgp;
308 : 1382592 : std::vector< tk::real > wgp;
309 : :
310 [ + - ]: 691296 : coordgp[0].resize( ng );
311 [ + - ]: 691296 : coordgp[1].resize( ng );
312 [ + - ]: 691296 : coordgp[2].resize( ng );
313 [ + - ]: 691296 : wgp.resize( ng );
314 : :
315 [ + - ]: 691296 : tk::GaussQuadratureTet( ng, coordgp, wgp );
316 : :
317 : : // Local degree of freedom
318 : 691296 : auto dof_el = ndofel[e];
319 : :
320 [ + - ]: 691296 : auto vole = geoElem(e, 0);
321 : :
322 : : // Loop over quadrature points in element e
323 [ + + ]: 1382592 : for (std::size_t igp=0; igp<ng; ++igp)
324 : : {
325 : : // Compute the basis function
326 [ + - ]: 1382592 : auto B = tk::eval_basis( dof_el, coordgp[0][igp], coordgp[1][igp],
327 : : coordgp[2][igp] );
328 : :
329 : 691296 : auto w = wgp[igp] * vole;
330 : :
331 [ + - ]: 1382592 : auto state = tk::eval_state( m_ncomp, rdof, dof_el, e, unk, B );
332 : :
333 : : // Mixture state at quadrature point
334 [ + - ]: 1382592 : Mixture mixgp(nspec, state, m_mat_blk);
335 : :
336 : : // Mixture density at quadrature point
337 : 691296 : tk::real rhob = mixgp.get_mix_density();
338 : :
339 : : // velocity vector at quadrature point
340 : : std::array< tk::real, 3 >
341 : 691296 : vel{ state[multispecies::momentumIdx(nspec, 0)]/rhob,
342 : 691296 : state[multispecies::momentumIdx(nspec, 1)]/rhob,
343 : 691296 : state[multispecies::momentumIdx(nspec, 2)]/rhob };
344 : :
345 [ + - ]: 1382592 : std::vector< tk::real > pri(m_nprim, 0.0);
346 : :
347 : : // Evaluate mixture temperature at quadrature point
348 : 691296 : auto rhoE0 = state[multispecies::energyIdx(nspec, 0)];
349 : 691296 : pri[multispecies::temperatureIdx(nspec,0)] =
350 [ + - ]: 691296 : mixgp.temperature(rhob, vel[0], vel[1], vel[2], rhoE0, m_mat_blk);
351 : : // TODO: consider clipping temperature here
352 : :
353 [ + + ]: 1382592 : for(std::size_t k = 0; k < m_nprim; k++)
354 : : {
355 : 691296 : auto mark = k * ndof;
356 [ + + ]: 1382592 : for(std::size_t idof = 0; idof < dof_el; idof++)
357 : 691296 : R[mark+idof] += w * pri[k] * B[idof];
358 : : }
359 : : }
360 : :
361 : : // Update the DG solution of primitive variables
362 [ + + ]: 1382592 : for(std::size_t k = 0; k < m_nprim; k++)
363 : : {
364 : 691296 : auto mark = k * ndof;
365 : 691296 : auto rmark = k * rdof;
366 [ + + ]: 1382592 : for(std::size_t idof = 0; idof < dof_el; idof++)
367 : : {
368 [ + - ]: 691296 : prim(e, rmark+idof) = R[mark+idof] / (mass_m[idof]*vole);
369 : : }
370 : : }
371 : : }
372 : 6688 : }
373 : :
374 : : //! Clean up the state of trace materials for this PDE system. No-op.
375 : : // //! \param[in] t Physical time
376 : : // //! \param[in] geoElem Element geometry array
377 : : // //! \param[in,out] unk Array of unknowns
378 : : // //! \param[in,out] prim Array of primitives
379 : : // //! \param[in] nielem Number of internal elements
380 : 6600 : void cleanTraceMaterial( tk::real /*t*/,
381 : : const tk::Fields& /*geoElem*/,
382 : : tk::Fields& /*unk*/,
383 : : tk::Fields& /*prim*/,
384 : 6600 : std::size_t /*nielem*/ ) const {}
385 : :
386 : : //! Reconstruct second-order solution from first-order
387 : : //! \param[in] geoElem Element geometry array
388 : : //! \param[in] fd Face connectivity and boundary conditions object
389 : : //! \param[in] esup Elements-surrounding-nodes connectivity
390 : : //! \param[in] inpoel Element-node connectivity
391 : : //! \param[in] coord Array of nodal coordinates
392 : : //! \param[in,out] U Solution vector at recent time step
393 : : //! \param[in,out] P Vector of primitives at recent time step
394 : : //! \param[in] pref Indicator for p-adaptive algorithm
395 : : //! \param[in] ndofel Vector of local number of degrees of freedome
396 : 6600 : void reconstruct( tk::real,
397 : : const tk::Fields&,
398 : : const tk::Fields& geoElem,
399 : : const inciter::FaceData& fd,
400 : : const std::map< std::size_t, std::vector< std::size_t > >&
401 : : esup,
402 : : const std::vector< std::size_t >& inpoel,
403 : : const tk::UnsMesh::Coords& coord,
404 : : tk::Fields& U,
405 : : tk::Fields& P,
406 : : const bool pref,
407 : : const std::vector< std::size_t >& ndofel ) const
408 : : {
409 : 6600 : const auto rdof = g_inputdeck.get< tag::rdof >();
410 : 6600 : const auto ndof = g_inputdeck.get< tag::ndof >();
411 : :
412 : 6600 : bool is_p0p1(false);
413 [ + - ][ + - ]: 6600 : if (rdof == 4 && ndof == 1)
414 : 6600 : is_p0p1 = true;
415 : :
416 : 6600 : const auto nelem = fd.Esuel().size()/4;
417 : :
418 [ - + ][ - - ]: 6600 : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
419 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
420 [ - + ][ - - ]: 6600 : Assert( P.nprop() == rdof*m_nprim, "Number of components in primitive "
[ - - ][ - - ]
[ - - ]
421 : : "vector must equal "+ std::to_string(rdof*m_nprim) );
422 : :
423 : : //----- reconstruction of conserved quantities -----
424 : : //--------------------------------------------------
425 : :
426 [ + + ]: 688800 : for (std::size_t e=0; e<nelem; ++e)
427 : : {
428 : 1364400 : std::vector< std::size_t > vars;
429 : : // check if element is marked as p0p1
430 [ - + ][ - - ]: 682200 : if ( (pref && ndofel[e] == 1) || is_p0p1 ) {
[ + - ][ + - ]
431 : : // 1. specify how many variables need to be reconstructed
432 [ + + ][ + - ]: 4434300 : for (std::size_t c=0; c<m_ncomp; ++c) vars.push_back(c);
433 : :
434 : : // 2. solve 3x3 least-squares system
435 : : // Reconstruct second-order dofs in Taylor space using nodal-stencils
436 [ + - ]: 682200 : tk::recoLeastSqExtStencil( rdof, e, esup, inpoel, geoElem, U, vars );
437 : :
438 : : // 3. transform reconstructed derivatives to Dubiner dofs
439 [ + - ]: 682200 : tk::transform_P0P1( rdof, e, inpoel, coord, U, vars );
440 : : }
441 : : }
442 : :
443 : : //----- reconstruction of primitive quantities -----
444 : : //--------------------------------------------------
445 : : // For multispecies, conserved and primitive quantities are reconstructed
446 : : // separately.
447 : :
448 [ + + ]: 688800 : for (std::size_t e=0; e<nelem; ++e)
449 : : {
450 : : // There are two conditions that requires the reconstruction of the
451 : : // primitive variables:
452 : : // 1. p-adaptive is triggered and P0P1 scheme is applied to specific
453 : : // elements
454 : : // 2. p-adaptive is not triggered and P0P1 scheme is applied to the
455 : : // whole computation domain
456 [ - + ][ - - ]: 682200 : if ((pref && ndofel[e] == 1) || (!pref && is_p0p1)) {
[ + - ][ + - ]
[ + - ]
457 : 1364400 : std::vector< std::size_t > vars;
458 [ + + ][ + - ]: 1364400 : for (std::size_t c=0; c<m_nprim; ++c) vars.push_back(c);
459 : :
460 : : // 1.
461 : : // Reconstruct second-order dofs in Taylor space using nodal-stencils
462 [ + - ]: 682200 : tk::recoLeastSqExtStencil( rdof, e, esup, inpoel, geoElem, P, vars );
463 : :
464 : : // 2.
465 [ + - ]: 682200 : tk::transform_P0P1(rdof, e, inpoel, coord, P, vars);
466 : : }
467 : : }
468 : 6600 : }
469 : :
470 : : //! Limit second-order solution, and primitive quantities separately
471 : : // //! \param[in] pref Indicator for p-adaptive algorithm
472 : : //! \param[in] geoFace Face geometry array
473 : : //! \param[in] geoElem Element geometry array
474 : : //! \param[in] fd Face connectivity and boundary conditions object
475 : : //! \param[in] esup Elements-surrounding-nodes connectivity
476 : : //! \param[in] inpoel Element-node connectivity
477 : : //! \param[in] coord Array of nodal coordinates
478 : : //! \param[in] ndofel Vector of local number of degrees of freedome
479 : : // //! \param[in] gid Local->global node id map
480 : : // //! \param[in] bid Local chare-boundary node ids (value) associated to
481 : : // //! global node ids (key)
482 : : // //! \param[in] uNodalExtrm Chare-boundary nodal extrema for conservative
483 : : // //! variables
484 : : // //! \param[in] pNodalExtrm Chare-boundary nodal extrema for primitive
485 : : // //! variables
486 : : // //! \param[in] mtInv Inverse of Taylor mass matrix
487 : : //! \param[in,out] U Solution vector at recent time step
488 : : //! \param[in,out] P Vector of primitives at recent time step
489 : : //! \param[in,out] shockmarker Vector of shock-marker values
490 : 6600 : void limit( [[maybe_unused]] tk::real,
491 : : const bool /*pref*/,
492 : : const tk::Fields& geoFace,
493 : : const tk::Fields& geoElem,
494 : : const inciter::FaceData& fd,
495 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
496 : : const std::vector< std::size_t >& inpoel,
497 : : const tk::UnsMesh::Coords& coord,
498 : : const std::vector< std::size_t >& ndofel,
499 : : const std::vector< std::size_t >& /*gid*/,
500 : : const std::unordered_map< std::size_t, std::size_t >& /*bid*/,
501 : : const std::vector< std::vector<tk::real> >& /*uNodalExtrm*/,
502 : : const std::vector< std::vector<tk::real> >& /*pNodalExtrm*/,
503 : : const std::vector< std::vector<tk::real> >& /*mtInv*/,
504 : : tk::Fields& U,
505 : : tk::Fields& P,
506 : : std::vector< std::size_t >& shockmarker ) const
507 : : {
508 : 6600 : const auto limiter = g_inputdeck.get< tag::limiter >();
509 : 6600 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
510 : 6600 : const auto rdof = g_inputdeck.get< tag::rdof >();
511 : 6600 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
512 : :
513 : : // limit vectors of conserved and primitive quantities
514 [ + - ][ + - ]: 6600 : if (limiter == ctr::LimiterType::VERTEXBASEDP1 && rdof == 4)
515 : : {
516 [ + - ]: 13200 : VertexBasedMultiSpecies_P1( esup, inpoel, ndofel, fd.Esuel().size()/4,
517 : 6600 : m_mat_blk, fd, geoFace, geoElem, coord, flux, solidx, U, P, nspec,
518 : : shockmarker );
519 : : }
520 [ - - ][ - - ]: 0 : else if (limiter == ctr::LimiterType::VERTEXBASEDP1 && rdof == 10)
521 : : {
522 [ - - ]: 0 : VertexBasedMultiSpecies_P2( esup, inpoel, ndofel, fd.Esuel().size()/4,
523 : 0 : m_mat_blk, fd, geoFace, geoElem, coord, flux, solidx, U, P, nspec,
524 : : shockmarker );
525 : : }
526 [ - - ]: 0 : else if (limiter != ctr::LimiterType::NOLIMITER)
527 : : {
528 [ - - ][ - - ]: 0 : Throw("Limiter type not configured for multispecies.");
[ - - ]
529 : : }
530 : 6600 : }
531 : :
532 : : //! Apply CPL to the conservative variable solution for this PDE system
533 : : //! \param[in] prim Array of primitive variables
534 : : //! \param[in] geoElem Element geometry array
535 : : //! \param[in] inpoel Element-node connectivity
536 : : //! \param[in] coord Array of nodal coordinates
537 : : //! \param[in,out] unk Array of conservative variables
538 : : //! \param[in] nielem Number of internal elements
539 : : //! \details This function applies CPL to obtain consistent dofs for
540 : : //! conservative quantities based on the limited primitive quantities.
541 : : //! See appendix of paper: Pandare et al. (2023). On the Design of Stable,
542 : : //! Consistent, and Conservative High-Order Methods for Multi-Material
543 : : //! Hydrodynamics. J Comp Phys, 112313.
544 : 6600 : void CPL( const tk::Fields& prim,
545 : : const tk::Fields& geoElem,
546 : : const std::vector< std::size_t >& inpoel,
547 : : const tk::UnsMesh::Coords& coord,
548 : : tk::Fields& unk,
549 : : std::size_t nielem ) const
550 : : {
551 : 6600 : [[maybe_unused]] const auto rdof = g_inputdeck.get< tag::rdof >();
552 : 6600 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
553 : :
554 [ - + ][ - - ]: 6600 : Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
555 : : "vector and primitive vector at recent time step incorrect" );
556 [ - + ][ - - ]: 6600 : Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
557 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
558 [ - + ][ - - ]: 6600 : Assert( prim.nprop() == rdof*m_nprim, "Number of components in vector of "
[ - - ][ - - ]
[ - - ]
559 : : "primitive quantities must equal "+ std::to_string(rdof*m_nprim) );
560 : :
561 : 6600 : correctLimConservMultiSpecies(nielem, m_mat_blk, nspec, inpoel,
562 : : coord, geoElem, prim, unk);
563 : 6600 : }
564 : :
565 : : //! Return cell-average deformation gradient tensor. No-op.
566 : 0 : std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
567 : : const tk::Fields&,
568 : : std::size_t ) const
569 : 0 : { return {}; }
570 : :
571 : : //! Reset the high order solution for p-adaptive scheme
572 : : //! \param[in] fd Face connectivity and boundary conditions object
573 : : //! \param[in,out] unk Solution vector at recent time step
574 : : //! \param[in,out] prim Primitive vector at recent time step
575 : : //! \param[in] ndofel Vector of local number of degrees of freedome
576 : : //! \details This function reset the high order coefficient for p-adaptive
577 : : //! solution polynomials. Unlike compflow class, the high order of fv
578 : : //! solution will not be reset since p0p1 is the base scheme for
579 : : //! multi-species p-adaptive DG method.
580 : 0 : void resetAdapSol( const inciter::FaceData& fd,
581 : : tk::Fields& unk,
582 : : tk::Fields& prim,
583 : : const std::vector< std::size_t >& ndofel ) const
584 : : {
585 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
586 : 0 : const auto ncomp = unk.nprop() / rdof;
587 : 0 : const auto nprim = prim.nprop() / rdof;
588 : :
589 [ - - ]: 0 : for(std::size_t e = 0; e < fd.Esuel().size()/4; e++)
590 : : {
591 [ - - ]: 0 : if(ndofel[e] < 10)
592 : : {
593 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c)
594 : : {
595 : 0 : auto mark = c*rdof;
596 : 0 : unk(e, mark+4) = 0.0;
597 : 0 : unk(e, mark+5) = 0.0;
598 : 0 : unk(e, mark+6) = 0.0;
599 : 0 : unk(e, mark+7) = 0.0;
600 : 0 : unk(e, mark+8) = 0.0;
601 : 0 : unk(e, mark+9) = 0.0;
602 : : }
603 [ - - ]: 0 : for (std::size_t c=0; c<nprim; ++c)
604 : : {
605 : 0 : auto mark = c*rdof;
606 : 0 : prim(e, mark+4) = 0.0;
607 : 0 : prim(e, mark+5) = 0.0;
608 : 0 : prim(e, mark+6) = 0.0;
609 : 0 : prim(e, mark+7) = 0.0;
610 : 0 : prim(e, mark+8) = 0.0;
611 : 0 : prim(e, mark+9) = 0.0;
612 : : }
613 : : }
614 : : }
615 : 0 : }
616 : :
617 : : //! Compute right hand side
618 : : //! \param[in] t Physical time
619 : : //! \param[in] pref Indicator for p-adaptive algorithm
620 : : //! \param[in] geoFace Face geometry array
621 : : //! \param[in] geoElem Element geometry array
622 : : //! \param[in] fd Face connectivity and boundary conditions object
623 : : //! \param[in] inpoel Element-node connectivity
624 : : //! \param[in] coord Array of nodal coordinates
625 : : //! \param[in] U Solution vector at recent time step
626 : : //! \param[in] P Primitive vector at recent time step
627 : : //! \param[in] ndofel Vector of local number of degrees of freedom
628 : : //! \param[in] dt Delta time
629 : : //! \param[in,out] R Right-hand side vector computed
630 : 6600 : void rhs( tk::real t,
631 : : const bool pref,
632 : : const tk::Fields& geoFace,
633 : : const tk::Fields& geoElem,
634 : : const inciter::FaceData& fd,
635 : : const std::vector< std::size_t >& inpoel,
636 : : const std::vector< std::unordered_set< std::size_t > >&,
637 : : const tk::UnsMesh::Coords& coord,
638 : : const tk::Fields& U,
639 : : const tk::Fields& P,
640 : : const std::vector< std::size_t >& ndofel,
641 : : const tk::real dt,
642 : : tk::Fields& R ) const
643 : : {
644 : 6600 : const auto ndof = g_inputdeck.get< tag::ndof >();
645 : 6600 : const auto rdof = g_inputdeck.get< tag::rdof >();
646 : 6600 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
647 : :
648 : 6600 : const auto nelem = fd.Esuel().size()/4;
649 : :
650 [ - + ][ - - ]: 6600 : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
651 : : "vector and primitive vector at recent time step incorrect" );
652 [ - + ][ - - ]: 6600 : Assert( U.nunk() == R.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
653 : : "vector and right-hand side at recent time step incorrect" );
654 [ - + ][ - - ]: 6600 : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
655 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
656 [ - + ][ - - ]: 6600 : Assert( P.nprop() == rdof*m_nprim, "Number of components in primitive "
[ - - ][ - - ]
[ - - ]
657 : : "vector must equal "+ std::to_string(rdof*m_nprim) );
658 [ - + ][ - - ]: 6600 : Assert( R.nprop() == ndof*m_ncomp, "Number of components in right-hand "
[ - - ][ - - ]
[ - - ]
659 : : "side vector must equal "+ std::to_string(ndof*m_ncomp) );
660 [ - + ][ - - ]: 6600 : Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
[ - - ][ - - ]
661 : : "Mismatch in inpofa size" );
662 : :
663 : : // set rhs to zero
664 [ + - ]: 6600 : R.fill(0.0);
665 : :
666 : : // empty vector for non-conservative terms. This vector is unused for
667 : : // multi-species flow since, there are no non-conservative terms
668 : : // in the system of PDEs.
669 : 13200 : std::vector< std::vector< tk::real > > riemannDeriv;
670 : :
671 : 13200 : std::vector< std::vector< tk::real > > vriem;
672 : 13200 : std::vector< std::vector< tk::real > > riemannLoc;
673 : :
674 : : // configure a no-op lambda for prescribed velocity
675 : 1752300 : auto velfn = []( ncomp_t, tk::real, tk::real, tk::real, tk::real ){
676 : 1752300 : return tk::VelFn::result_type(); };
677 : :
678 : : // compute internal surface flux integrals
679 [ + - ]: 13200 : tk::surfInt( pref, 1, m_mat_blk, t, ndof, rdof, inpoel, solidx,
680 [ + - ]: 6600 : coord, fd, geoFace, geoElem, m_riemann, velfn, U, P, ndofel,
681 : : dt, R, riemannDeriv );
682 : :
683 : : // compute optional source term
684 [ + - ][ + - ]: 6600 : tk::srcInt( m_mat_blk, t, ndof, fd.Esuel().size()/4, inpoel,
685 : : coord, geoElem, Problem::src, ndofel, R );
686 : :
687 [ - + ]: 6600 : if(ndof > 1)
688 : : // compute volume integrals
689 [ - - ][ - - ]: 0 : tk::volInt( 1, t, m_mat_blk, ndof, rdof, nelem, inpoel, coord, geoElem,
[ - - ]
690 : : flux, velfn, U, P, ndofel, R );
691 : :
692 : : // compute boundary surface flux integrals
693 [ + + ]: 52800 : for (const auto& b : m_bc)
694 [ + - ][ + - ]: 92400 : tk::bndSurfInt( pref, 1, m_mat_blk, ndof, rdof, std::get<0>(b), fd,
695 : 46200 : geoFace, geoElem, inpoel, coord, t, m_riemann, velfn,
696 : 46200 : std::get<1>(b), U, P, ndofel, R, riemannDeriv );
697 : :
698 : : // compute external (energy) sources
699 : : //m_physics.physSrc(nspec, t, geoElem, {}, R, {});
700 : 6600 : }
701 : :
702 : : //! Evaluate the adaptive indicator and mark the ndof for each element
703 : : //! \param[in] nunk Number of unknowns
704 : : //! \param[in] coord Array of nodal coordinates
705 : : //! \param[in] inpoel Element-node connectivity
706 : : //! \param[in] fd Face connectivity and boundary conditions object
707 : : //! \param[in] unk Array of unknowns
708 : : // //! \param[in] prim Array of primitive quantities
709 : : //! \param[in] indicator p-refinement indicator type
710 : : //! \param[in] ndof Number of degrees of freedom in the solution
711 : : //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
712 : : //! \param[in] tolref Tolerance for p-refinement
713 : : //! \param[in,out] ndofel Vector of local number of degrees of freedome
714 : 0 : void eval_ndof( std::size_t nunk,
715 : : [[maybe_unused]] const tk::UnsMesh::Coords& coord,
716 : : [[maybe_unused]] const std::vector< std::size_t >& inpoel,
717 : : const inciter::FaceData& fd,
718 : : const tk::Fields& unk,
719 : : const tk::Fields& /*prim*/,
720 : : inciter::ctr::PrefIndicatorType indicator,
721 : : std::size_t ndof,
722 : : std::size_t ndofmax,
723 : : tk::real tolref,
724 : : std::vector< std::size_t >& ndofel ) const
725 : : {
726 : 0 : const auto& esuel = fd.Esuel();
727 : :
728 [ - - ]: 0 : if(indicator == inciter::ctr::PrefIndicatorType::SPECTRAL_DECAY)
729 : 0 : spectral_decay(1, nunk, esuel, unk, ndof, ndofmax, tolref, ndofel);
730 : : else
731 [ - - ][ - - ]: 0 : Throw( "No such adaptive indicator type" );
[ - - ]
732 : 0 : }
733 : :
734 : : //! Compute the minimum time step size
735 : : //! \param[in] fd Face connectivity and boundary conditions object
736 : : //! \param[in] geoFace Face geometry array
737 : : //! \param[in] geoElem Element geometry array
738 : : // //! \param[in] ndofel Vector of local number of degrees of freedom
739 : : //! \param[in] U Solution vector at recent time step
740 : : //! \param[in] P Vector of primitive quantities at recent time step
741 : : //! \param[in] nielem Number of internal elements
742 : : //! \return Minimum time step size
743 : : //! \details The allowable dt is calculated by looking at the maximum
744 : : //! wave-speed in elements surrounding each face, times the area of that
745 : : //! face. Once the maximum of this quantity over the mesh is determined,
746 : : //! the volume of each cell is divided by this quantity. A minimum of this
747 : : //! ratio is found over the entire mesh, which gives the allowable dt.
748 : 0 : tk::real dt( const std::array< std::vector< tk::real >, 3 >&,
749 : : const std::vector< std::size_t >&,
750 : : const inciter::FaceData& fd,
751 : : const tk::Fields& geoFace,
752 : : const tk::Fields& geoElem,
753 : : const std::vector< std::size_t >& /*ndofel*/,
754 : : const tk::Fields& U,
755 : : const tk::Fields& P,
756 : : const std::size_t nielem ) const
757 : : {
758 : 0 : const auto ndof = g_inputdeck.get< tag::ndof >();
759 : 0 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
760 : :
761 : 0 : auto mindt = timeStepSizeMultiSpecies( m_mat_blk, fd.Esuf(), geoFace,
762 : : geoElem, nielem, nspec, U, P);
763 : :
764 : : //if (viscous)
765 : : // mindt = std::min(mindt, timeStepSizeViscousFV(geoElem, nielem, nspec, U));
766 : : //mindt = std::min(mindt, m_physics.dtRestriction(geoElem, nielem, {}));
767 : :
768 : 0 : tk::real dgp = 0.0;
769 [ - - ]: 0 : if (ndof == 4)
770 : : {
771 : 0 : dgp = 1.0;
772 : : }
773 [ - - ]: 0 : else if (ndof == 10)
774 : : {
775 : 0 : dgp = 2.0;
776 : : }
777 : :
778 : : // Scale smallest dt with CFL coefficient and the CFL is scaled by (2*p+1)
779 : : // where p is the order of the DG polynomial by linear stability theory.
780 : 0 : mindt /= (2.0*dgp + 1.0);
781 : 0 : return mindt;
782 : : }
783 : :
784 : : //! Balances elastic energy after plastic update. Not implemented here.
785 : : // //! \param[in] e Element number
786 : : // //! \param[in] x_star Stiff variables before implicit update
787 : : // //! \param[in] x Stiff variables after implicit update
788 : : // //! \param[in] U Field of conserved variables
789 : 0 : void balance_plastic_energy( std::size_t /*e*/,
790 : : std::vector< tk::real > /*x_star*/,
791 : : std::vector< tk::real > /*x*/,
792 : 0 : tk::Fields& /*U*/ ) const {}
793 : :
794 : : //! Compute stiff terms for a single element. No-op until chem sources added
795 : : // //! \param[in] e Element number
796 : : // //! \param[in] geoElem Element geometry array
797 : : // //! \param[in] inpoel Element-node connectivity
798 : : // //! \param[in] coord Array of nodal coordinates
799 : : // //! \param[in] U Solution vector at recent time step
800 : : // //! \param[in] P Primitive vector at recent time step
801 : : // //! \param[in] ndofel Vector of local number of degrees of freedom
802 : : // //! \param[in,out] R Right-hand side vector computed
803 : 0 : void stiff_rhs( std::size_t /*e*/,
804 : : const tk::Fields& /*geoElem*/,
805 : : const std::vector< std::size_t >& /*inpoel*/,
806 : : const tk::UnsMesh::Coords& /*coord*/,
807 : : const tk::Fields& /*U*/,
808 : : const tk::Fields& /*P*/,
809 : : const std::vector< std::size_t >& /*ndofel*/,
810 : 0 : tk::Fields& /*R*/ ) const {}
811 : :
812 : : //! Extract the velocity field at cell nodes. Currently unused.
813 : : // //! \param[in] U Solution vector at recent time step
814 : : // //! \param[in] N Element node indices
815 : : //! \return Array of the four values of the velocity field
816 : : std::array< std::array< tk::real, 4 >, 3 >
817 : : velocity( const tk::Fields& /*U*/,
818 : : const std::array< std::vector< tk::real >, 3 >&,
819 : : const std::array< std::size_t, 4 >& /*N*/ ) const
820 : : {
821 : : std::array< std::array< tk::real, 4 >, 3 > v;
822 : : return v;
823 : : }
824 : :
825 : : //! Return a map that associates user-specified strings to functions
826 : : //! \return Map that associates user-specified strings to functions that
827 : : //! compute relevant quantities to be output to file
828 : 352 : std::map< std::string, tk::GetVarFn > OutVarFn() const
829 : 352 : { return MultiSpeciesOutVarFn(); }
830 : :
831 : : //! Return analytic field names to be output to file
832 : : //! \return Vector of strings labelling analytic fields output in file
833 : 0 : std::vector< std::string > analyticFieldNames() const {
834 : 0 : auto nspec = g_inputdeck.get< eq, tag::nspec >();
835 : :
836 : 0 : return MultiSpeciesFieldNames(nspec);
837 : : }
838 : :
839 : : //! Return time history field names to be output to file
840 : : //! \return Vector of strings labelling time history fields output in file
841 : 0 : std::vector< std::string > histNames() const {
842 : 0 : return MultiSpeciesHistNames();
843 : : }
844 : :
845 : : //! Return surface field names to be output to file
846 : : //! \return Vector of strings labelling surface fields output in file
847 : 176 : std::vector< std::string > surfNames() const
848 : : {
849 : 176 : return MultiSpeciesSurfNames();
850 : : }
851 : :
852 : : //! Return surface field output going to file
853 : : std::vector< std::vector< tk::real > >
854 : 176 : surfOutput( const inciter::FaceData& fd,
855 : : const tk::Fields& U,
856 : : const tk::Fields& P ) const
857 : : {
858 : 176 : const auto rdof = g_inputdeck.get< tag::rdof >();
859 : 176 : const auto nspec = g_inputdeck.get< eq, tag::nspec >();
860 : :
861 : 176 : return MultiSpeciesSurfOutput( nspec, rdof, fd, U, P );
862 : : }
863 : :
864 : : //! Return time history field output evaluated at time history points
865 : : //! \param[in] h History point data
866 : : //! \param[in] inpoel Element-node connectivity
867 : : //! \param[in] coord Array of nodal coordinates
868 : : //! \param[in] U Array of unknowns
869 : : //! \param[in] P Array of primitive quantities
870 : : //! \return Vector of time history output of bulk flow quantities (density,
871 : : //! velocity, total energy, and pressure) evaluated at time history points
872 : : std::vector< std::vector< tk::real > >
873 : 0 : histOutput( const std::vector< HistData >& h,
874 : : const std::vector< std::size_t >& inpoel,
875 : : const tk::UnsMesh::Coords& coord,
876 : : const tk::Fields& U,
877 : : const tk::Fields& P ) const
878 : : {
879 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
880 : 0 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
881 : :
882 : 0 : const auto& x = coord[0];
883 : 0 : const auto& y = coord[1];
884 : 0 : const auto& z = coord[2];
885 : :
886 [ - - ]: 0 : std::vector< std::vector< tk::real > > Up(h.size());
887 : :
888 : 0 : std::size_t j = 0;
889 [ - - ]: 0 : for (const auto& p : h) {
890 : 0 : auto e = p.get< tag::elem >();
891 : 0 : auto chp = p.get< tag::coord >();
892 : :
893 : : // Evaluate inverse Jacobian
894 : 0 : std::array< std::array< tk::real, 3>, 4 > cp{{
895 : : {{ x[inpoel[4*e ]], y[inpoel[4*e ]], z[inpoel[4*e ]] }},
896 : 0 : {{ x[inpoel[4*e+1]], y[inpoel[4*e+1]], z[inpoel[4*e+1]] }},
897 : 0 : {{ x[inpoel[4*e+2]], y[inpoel[4*e+2]], z[inpoel[4*e+2]] }},
898 : 0 : {{ x[inpoel[4*e+3]], y[inpoel[4*e+3]], z[inpoel[4*e+3]] }} }};
899 : 0 : auto J = tk::inverseJacobian( cp[0], cp[1], cp[2], cp[3] );
900 : :
901 : : // evaluate solution at history-point
902 : 0 : std::array< tk::real, 3 > dc{{chp[0]-cp[0][0], chp[1]-cp[0][1],
903 : 0 : chp[2]-cp[0][2]}};
904 [ - - ]: 0 : auto B = tk::eval_basis(rdof, tk::dot(J[0],dc), tk::dot(J[1],dc),
905 : 0 : tk::dot(J[2],dc));
906 [ - - ]: 0 : auto uhp = eval_state(m_ncomp, rdof, rdof, e, U, B);
907 [ - - ]: 0 : auto php = eval_state(m_nprim, rdof, rdof, e, P, B);
908 : :
909 : : // Mixture calculations, initialized
910 [ - - ]: 0 : Mixture mix(nspec, uhp, m_mat_blk);
911 : :
912 : : // store solution in history output vector
913 [ - - ]: 0 : Up[j].resize(6+nspec, 0.0);
914 : 0 : Up[j][0] = mix.get_mix_density();
915 : 0 : Up[j][1] = uhp[multispecies::momentumIdx(nspec,0)]/Up[j][0];
916 : 0 : Up[j][2] = uhp[multispecies::momentumIdx(nspec,1)]/Up[j][0];
917 : 0 : Up[j][3] = uhp[multispecies::momentumIdx(nspec,2)]/Up[j][0];
918 : 0 : Up[j][4] = uhp[multispecies::energyIdx(nspec,0)];
919 [ - - ]: 0 : Up[j][5] = mix.pressure( Up[j][0],
920 : : php[multispecies::temperatureIdx(nspec,0)] );
921 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k) {
922 : 0 : Up[j][6+k] = uhp[multispecies::densityIdx(nspec,k)]/Up[j][0];
923 : : }
924 : 0 : ++j;
925 : : }
926 : :
927 : 0 : return Up;
928 : : }
929 : :
930 : : //! Return names of integral variables to be output to diagnostics file
931 : : //! \return Vector of strings labelling integral variables output
932 : 6 : std::vector< std::string > names() const
933 : : {
934 : 6 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
935 : 6 : return MultiSpeciesDiagNames(nspec);
936 : : }
937 : :
938 : : //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
939 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
940 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
941 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
942 : : //! \param[in] t Physical time at which to evaluate the analytic solution
943 : : //! \return Vector of analytic solution at given location and time
944 : : std::vector< tk::real >
945 : 0 : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
946 : 0 : { return Problem::analyticSolution( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
947 : :
948 : : //! Return analytic solution for conserved variables
949 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
950 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
951 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
952 : : //! \param[in] t Physical time at which to evaluate the analytic solution
953 : : //! \return Vector of analytic solution at given location and time
954 : : std::vector< tk::real >
955 : 227400 : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
956 : 227400 : { return Problem::initialize( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
957 : :
958 : : //! Return cell-averaged specific total energy for an element
959 : : //! \param[in] e Element id for which total energy is required
960 : : //! \param[in] unk Vector of conserved quantities
961 : : //! \return Cell-averaged specific total energy for given element
962 : 227400 : tk::real sp_totalenergy(std::size_t e, const tk::Fields& unk) const
963 : : {
964 : 227400 : const auto rdof = g_inputdeck.get< tag::rdof >();
965 : 227400 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
966 : :
967 : 227400 : return unk(e, multispecies::energyDofIdx(nspec,0,rdof,0));
968 : : }
969 : :
970 : : private:
971 : : //! Physics policy
972 : : const Physics m_physics;
973 : : //! Number of components in this PDE system
974 : : const ncomp_t m_ncomp;
975 : : //! Number of primitive quantities stored in this PDE system
976 : : const ncomp_t m_nprim;
977 : : //! Riemann solver
978 : : tk::RiemannFluxFn m_riemann;
979 : : //! BC configuration
980 : : BCStateFn m_bc;
981 : : //! EOS material block
982 : : std::vector< EOS > m_mat_blk;
983 : :
984 : : //! Evaluate conservative part of physical flux function for this PDE system
985 : : //! \param[in] ncomp Number of scalar components in this PDE system
986 : : //! \param[in] ugp Numerical solution at the Gauss point at which to
987 : : //! evaluate the flux
988 : : //! \return Flux vectors for all components in this PDE system
989 : : //! \note The function signature must follow tk::FluxFn
990 : : static tk::FluxFn::result_type
991 : 0 : flux( [[maybe_unused]] ncomp_t ncomp,
992 : : const std::vector< EOS >& mat_blk,
993 : : const std::vector< tk::real >& ugp,
994 : : const std::vector< std::array< tk::real, 3 > >& )
995 : : {
996 : 0 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
997 : :
998 [ - - ]: 0 : std::vector< std::array< tk::real, 3 > > fl( ugp.size() );
999 : :
1000 [ - - ]: 0 : Mixture mix(nspec, ugp, mat_blk);
1001 : 0 : auto rhob = mix.get_mix_density();
1002 : :
1003 : 0 : std::array< tk::real, 3 > u{{
1004 : 0 : ugp[multispecies::momentumIdx(nspec,0)] / rhob,
1005 : 0 : ugp[multispecies::momentumIdx(nspec,1)] / rhob,
1006 : 0 : ugp[multispecies::momentumIdx(nspec,2)] / rhob }};
1007 [ - - ]: 0 : auto p = mix.pressure(rhob,
1008 : 0 : ugp[ncomp+multispecies::temperatureIdx(nspec,0)]);
1009 : :
1010 : : // density flux
1011 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k) {
1012 : 0 : auto idx = multispecies::densityIdx(nspec, k);
1013 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
1014 : 0 : fl[idx][j] = ugp[idx] * u[j];
1015 : : }
1016 : : }
1017 : :
1018 : : // momentum flux
1019 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
1020 : 0 : auto idx = multispecies::momentumIdx(nspec,i);
1021 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
1022 : 0 : fl[idx][j] = ugp[idx] * u[j];
1023 [ - - ]: 0 : if (i == j) fl[idx][j] += p;
1024 : : }
1025 : : }
1026 : :
1027 : : // energy flux
1028 : 0 : auto idx = multispecies::energyIdx(nspec,0);
1029 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
1030 : 0 : fl[idx][j] = u[j] * (ugp[idx] + p);
1031 : : }
1032 : :
1033 : 0 : return fl;
1034 : : }
1035 : :
1036 : : //! \brief Boundary state function providing the left and right state of a
1037 : : //! face at Dirichlet boundaries
1038 : : //! \param[in] ncomp Number of scalar components in this PDE system
1039 : : //! \param[in] mat_blk EOS material block
1040 : : //! \param[in] ul Left (domain-internal) state
1041 : : //! \param[in] x X-coordinate at which to compute the states
1042 : : //! \param[in] y Y-coordinate at which to compute the states
1043 : : //! \param[in] z Z-coordinate at which to compute the states
1044 : : //! \param[in] t Physical time
1045 : : //! \return Left and right states for all scalar components in this PDE
1046 : : //! system
1047 : : //! \note The function signature must follow tk::StateFn.
1048 : : static tk::StateFn::result_type
1049 : 2250 : dirichlet( ncomp_t ncomp,
1050 : : const std::vector< EOS >& mat_blk,
1051 : : const std::vector< tk::real >& ul, tk::real x, tk::real y,
1052 : : tk::real z, tk::real t, const std::array< tk::real, 3 >& )
1053 : : {
1054 : 2250 : return {{ ul, Problem::initialize( ncomp, mat_blk, x, y, z, t ) }};
1055 : : }
1056 : :
1057 : : // Other boundary condition types that do not depend on "Problem" should be
1058 : : // added in BCFunctions.hpp
1059 : : };
1060 : :
1061 : : } // dg::
1062 : :
1063 : : } // inciter::
1064 : :
1065 : : #endif // DGMultiSpecies_h
|
