Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/MultiMat/DGMultiMat.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-material flow using discontinuous Galerkin
9 : : finite elements
10 : : \details This file implements calls to the physics operators governing
11 : : compressible multi-material flow using discontinuous Galerkin
12 : : discretizations.
13 : : */
14 : : // *****************************************************************************
15 : : #ifndef DGMultiMat_h
16 : : #define DGMultiMat_h
17 : :
18 : : #include <cmath>
19 : : #include <algorithm>
20 : : #include <unordered_set>
21 : : #include <map>
22 : : #include <array>
23 : :
24 : : #include "Macro.hpp"
25 : : #include "Exception.hpp"
26 : : #include "Vector.hpp"
27 : : #include "ContainerUtil.hpp"
28 : : #include "UnsMesh.hpp"
29 : : #include "Inciter/InputDeck/InputDeck.hpp"
30 : : #include "Integrate/Basis.hpp"
31 : : #include "Integrate/Quadrature.hpp"
32 : : #include "Integrate/Initialize.hpp"
33 : : #include "Integrate/Mass.hpp"
34 : : #include "Integrate/Surface.hpp"
35 : : #include "Integrate/Boundary.hpp"
36 : : #include "Integrate/Volume.hpp"
37 : : #include "Integrate/MultiMatTerms.hpp"
38 : : #include "Integrate/Source.hpp"
39 : : #include "RiemannFactory.hpp"
40 : : #include "EoS/EoS.hpp"
41 : : #include "MultiMat/MultiMatIndexing.hpp"
42 : : #include "Reconstruction.hpp"
43 : : #include "Limiter.hpp"
44 : : #include "Problem/FieldOutput.hpp"
45 : : #include "Problem/BoxInitialization.hpp"
46 : : #include "PrefIndicator.hpp"
47 : :
48 : : namespace inciter {
49 : :
50 : : extern ctr::InputDeck g_inputdeck;
51 : :
52 : : namespace dg {
53 : :
54 : : //! \brief MultiMat used polymorphically with tk::DGPDE
55 : : //! \details The template arguments specify policies and are used to configure
56 : : //! the behavior of the class. The policies are:
57 : : //! - Physics - physics configuration, see PDE/MultiMat/Physics.h
58 : : //! - Problem - problem configuration, see PDE/MultiMat/Problem.h
59 : : //! \note The default physics is velocity equilibrium (veleq), set in
60 : : //! inciter::deck::check_multimat()
61 : : template< class Physics, class Problem >
62 [ - - ][ - - ]: 30 : class MultiMat {
[ + - ][ + - ]
[ + - ][ - - ]
63 : :
64 : : private:
65 : : using eq = tag::multimat;
66 : :
67 : : public:
68 : : //! Constructor
69 : : //! \param[in] c Equation system index (among multiple systems configured)
70 : 30 : explicit MultiMat( ncomp_t c ) :
71 : : m_system( c ),
72 : : m_ncomp( g_inputdeck.get< tag::component, eq >().at(c) ),
73 : 30 : m_offset( g_inputdeck.get< tag::component >().offset< eq >(c) ),
74 : : m_riemann( tk::cref_find( multimatRiemannSolvers(),
75 [ - + ][ - + ]: 90 : g_inputdeck.get< tag::param, tag::multimat, tag::flux >().at(m_system) ) )
[ + - ][ + - ]
[ + - ]
76 : : {
77 : : // associate boundary condition configurations with state functions
78 [ + - ][ + - ]: 210 : brigand::for_each< ctr::bc::Keys >( ConfigBC< eq >( m_system, m_bc,
[ + + ][ + - ]
[ - - ][ - - ]
79 : : { dirichlet
80 : : , symmetry
81 : : , invalidBC // Inlet BC not implemented
82 : : , invalidBC // Outlet BC not implemented
83 : : , subsonicOutlet
84 : : , extrapolate } ) );
85 : 30 : }
86 : :
87 : : //! Find the number of primitive quantities required for this PDE system
88 : : //! \return The number of primitive quantities required to be stored for
89 : : //! this PDE system
90 : : std::size_t nprim() const
91 : : {
92 : 47868523 : const auto nmat =
93 : : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[m_system];
94 : : // multimat needs individual material pressures and velocities currently
95 [ - - ][ - - ]: 43314481 : return (nmat+3);
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ + - ]
[ + - ][ + - ]
[ + + ][ + + ]
[ - - ][ - - ]
[ - - ][ - - ]
[ + + ][ + + ]
[ - - ][ + - ]
[ + - ][ + - ]
[ + + ][ + + ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
96 : : }
97 : :
98 : : //! Find the number of materials set up for this PDE system
99 : : //! \return The number of materials set up for this PDE system
100 : : std::size_t nmat() const
101 : : {
102 : 0 : const auto nmat =
103 : : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[m_system];
104 : : return nmat;
105 : : }
106 : :
107 : : //! Assign number of DOFs per equation in the PDE system
108 : : //! \param[in,out] numEqDof Array storing number of Dofs for each PDE
109 : : //! equation
110 : 65 : void numEquationDofs(std::vector< std::size_t >& numEqDof) const
111 : : {
112 : : // all equation-dofs initialized to ndofs first
113 [ + + ]: 665 : for (std::size_t i=0; i<m_ncomp; ++i) {
114 : 600 : numEqDof.push_back(g_inputdeck.get< tag::discr, tag::ndof >());
115 : : }
116 : :
117 : : // volume fractions are P0Pm (ndof = 1) if interface reconstruction is used
118 : 65 : const auto nmat =
119 : 65 : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[m_system];
120 [ + + ]: 200 : for (std::size_t k=0; k<nmat; ++k) {
121 [ + + ]: 135 : if (g_inputdeck.get< tag::param, tag::multimat, tag::intsharp >
122 [ + + ]: 135 : ()[m_system] > 0) numEqDof[volfracIdx(nmat, k)] = 1;
123 : : }
124 : 65 : }
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 : : void IcBoxElems( const tk::Fields& geoElem,
132 : : std::size_t nielem,
133 : : std::vector< std::unordered_set< std::size_t > >& inbox ) const
134 : : {
135 : 65 : tk::BoxElems< eq >(m_system, geoElem, nielem, inbox);
136 : : }
137 : :
138 : : //! Initalize the compressible flow equations, prepare for time integration
139 : : //! \param[in] L Block diagonal mass matrix
140 : : //! \param[in] inpoel Element-node connectivity
141 : : //! \param[in] coord Array of nodal coordinates
142 : : //! \param[in] inbox List of elements at which box user ICs are set for
143 : : //! each IC box
144 : : //! \param[in,out] unk Array of unknowns
145 : : //! \param[in] t Physical time
146 : : //! \param[in] nielem Number of internal elements
147 [ + - ]: 65 : void initialize( const tk::Fields& L,
148 : : const std::vector< std::size_t >& inpoel,
149 : : const tk::UnsMesh::Coords& coord,
150 : : const std::vector< std::unordered_set< std::size_t > >& inbox,
151 : : tk::Fields& unk,
152 : : tk::real t,
153 : : const std::size_t nielem ) const
154 : : {
155 [ + - ]: 65 : tk::initialize( m_system, m_ncomp, m_offset, L, inpoel, coord,
156 : : Problem::initialize, unk, t, nielem );
157 : :
158 : 65 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
159 : : const auto& ic = g_inputdeck.get< tag::param, eq, tag::ic >();
160 : : const auto& icbox = ic.get< tag::box >();
161 : :
162 : : // Set initial conditions inside user-defined IC box
163 : 65 : std::vector< tk::real > s(m_ncomp, 0.0);
164 [ + + ]: 22511 : for (std::size_t e=0; e<nielem; ++e) {
165 [ + - ]: 22446 : if (icbox.size() > m_system) {
166 : : std::size_t bcnt = 0;
167 [ - + ]: 22446 : for (const auto& b : icbox[m_system]) { // for all boxes
168 [ - - ][ - - ]: 0 : if (inbox.size() > bcnt && inbox[bcnt].find(e) != inbox[bcnt].end())
169 : : {
170 [ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) {
171 : 0 : auto mark = c*rdof;
172 : 0 : s[c] = unk(e,mark,m_offset);
173 : : // set high-order DOFs to zero
174 [ - - ]: 0 : for (std::size_t i=1; i<rdof; ++i)
175 : 0 : unk(e,mark+i,m_offset) = 0.0;
176 : : }
177 [ - - ]: 0 : initializeBox( m_system, 1.0, t, b, s );
178 : : // store box-initialization in solution vector
179 [ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) {
180 : 0 : auto mark = c*rdof;
181 : 0 : unk(e,mark,m_offset) = s[c];
182 : : }
183 : : }
184 : 0 : ++bcnt;
185 : : }
186 : : }
187 : : }
188 : 65 : }
189 : :
190 : : //! Compute the left hand side block-diagonal mass matrix
191 : : //! \param[in] geoElem Element geometry array
192 : : //! \param[in,out] l Block diagonal mass matrix
193 : : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const {
194 : 65 : const auto ndof = g_inputdeck.get< tag::discr, tag::ndof >();
195 : 65 : tk::mass( m_ncomp, m_offset, ndof, geoElem, l );
196 : : }
197 : :
198 : : //! Update the interface cells to first order dofs
199 : : //! \param[in] unk Array of unknowns
200 : : //! \param[in] nielem Number of internal elements
201 : : // //! \param[in,out] ndofel Array of dofs
202 : : //! \details This function resets the high-order terms in interface cells.
203 : 5250 : void updateInterfaceCells( tk::Fields& unk,
204 : : std::size_t nielem,
205 : : std::vector< std::size_t >& /*ndofel*/ ) const
206 : : {
207 : 5250 : auto intsharp =
208 [ + + ]: 5250 : g_inputdeck.get< tag::param, tag::multimat, tag::intsharp >()[m_system];
209 : : // If this cell is not material interface, return this function
210 [ + + ]: 5250 : if(not intsharp) return;
211 : :
212 : 375 : auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
213 : 375 : auto nmat =
214 : : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[m_system];
215 : :
216 [ + + ]: 227775 : for (std::size_t e=0; e<nielem; ++e) {
217 : 227400 : std::vector< std::size_t > matInt(nmat, 0);
218 [ + - ][ - - ]: 227400 : std::vector< tk::real > alAvg(nmat, 0.0);
219 [ + + ]: 682200 : for (std::size_t k=0; k<nmat; ++k)
220 : 454800 : alAvg[k] = unk(e, volfracDofIdx(nmat,k,rdof,0), m_offset);
221 [ + - ]: 227400 : auto intInd = interfaceIndicator(nmat, alAvg, matInt);
222 : :
223 : : // interface cells cannot be high-order
224 [ + + ]: 227400 : if (intInd) {
225 : : //ndofel[e] = 1;
226 [ + + ]: 32478 : for (std::size_t k=0; k<nmat; ++k) {
227 [ + - ]: 21652 : if (matInt[k]) {
228 [ + + ]: 86608 : for (std::size_t i=1; i<rdof; ++i) {
229 : 64956 : unk(e, densityDofIdx(nmat,k,rdof,i), m_offset) = 0.0;
230 : 64956 : unk(e, energyDofIdx(nmat,k,rdof,i), m_offset) = 0.0;
231 : : }
232 : : }
233 : : }
234 [ + + ]: 43304 : for (std::size_t idir=0; idir<3; ++idir) {
235 [ + + ]: 129912 : for (std::size_t i=1; i<rdof; ++i) {
236 : 97434 : unk(e, momentumDofIdx(nmat,idir,rdof,i), m_offset) = 0.0;
237 : : }
238 : : }
239 : : }
240 : : }
241 : : }
242 : :
243 : : //! Update the primitives for this PDE system
244 : : //! \param[in] unk Array of unknowns
245 : : //! \param[in] L The left hand side block-diagonal mass matrix
246 : : //! \param[in] geoElem Element geometry array
247 : : //! \param[in,out] prim Array of primitives
248 : : //! \param[in] nielem Number of internal elements
249 : : //! \details This function computes and stores the dofs for primitive
250 : : //! quantities, which are required for obtaining reconstructed states used
251 : : //! in the Riemann solver. See /PDE/Riemann/AUSM.hpp, where the
252 : : //! normal velocity for advection is calculated from independently
253 : : //! reconstructed velocities.
254 : 5315 : void updatePrimitives( const tk::Fields& unk,
255 : : const tk::Fields& L,
256 : : const tk::Fields& geoElem,
257 : : tk::Fields& prim,
258 : : std::size_t nielem ) const
259 : : {
260 : 5315 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
261 : 5315 : const auto ndof = g_inputdeck.get< tag::discr, tag::ndof >();
262 : 5315 : const auto nmat =
263 : 5315 : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[m_system];
264 : :
265 : : Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
266 : : "vector and primitive vector at recent time step incorrect" );
267 : : Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
268 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
269 : : Assert( prim.nprop() == rdof*nprim(), "Number of components in vector of "
270 : : "primitive quantities must equal "+ std::to_string(rdof*nprim()) );
271 : : Assert( (g_inputdeck.get< tag::discr, tag::ndof >()) <= 4, "High-order "
272 : : "discretizations not set up for multimat updatePrimitives()" );
273 : :
274 [ + + ]: 2257661 : for (std::size_t e=0; e<nielem; ++e)
275 : : {
276 : 2252346 : std::vector< tk::real > R(nprim()*ndof, 0.0);
277 : :
278 [ + - ]: 2252346 : auto ng = tk::NGvol(ndof);
279 : :
280 : : // arrays for quadrature points
281 : : std::array< std::vector< tk::real >, 3 > coordgp;
282 : : std::vector< tk::real > wgp;
283 : :
284 [ + - ]: 2252346 : coordgp[0].resize( ng );
285 [ + - ]: 2252346 : coordgp[1].resize( ng );
286 [ + - ]: 2252346 : coordgp[2].resize( ng );
287 [ + - ]: 2252346 : wgp.resize( ng );
288 : :
289 [ + - ]: 2252346 : tk::GaussQuadratureTet( ng, coordgp, wgp );
290 : :
291 : : // Loop over quadrature points in element e
292 [ + + ]: 6348148 : for (std::size_t igp=0; igp<ng; ++igp)
293 : : {
294 : : // Compute the basis function
295 [ + - ]: 4095802 : auto B =
296 [ + - ]: 4095802 : tk::eval_basis( ndof, coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
297 : :
298 [ + - ]: 4095802 : auto w = wgp[igp] * geoElem(e, 0, 0);
299 : :
300 [ + - ][ - - ]: 4095802 : auto state = tk::eval_state( m_ncomp, 0, rdof, ndof, e, unk, B, {0, m_ncomp-1} );
301 : :
302 : : // bulk density at quadrature point
303 : : tk::real rhob(0.0);
304 [ + + ]: 13387592 : for (std::size_t k=0; k<nmat; ++k)
305 : 9291790 : rhob += state[densityIdx(nmat, k)];
306 : :
307 : : // velocity vector at quadrature point
308 : : std::array< tk::real, 3 >
309 [ + - ]: 4095802 : vel{ state[momentumIdx(nmat, 0)]/rhob,
310 : 4095802 : state[momentumIdx(nmat, 1)]/rhob,
311 : 4095802 : state[momentumIdx(nmat, 2)]/rhob };
312 : :
313 [ + - ][ - - ]: 4095802 : std::vector< tk::real > pri(nprim(), 0.0);
314 : :
315 : : // Evaluate material pressure at quadrature point
316 [ + + ]: 13387592 : for(std::size_t imat = 0; imat < nmat; imat++)
317 : : {
318 [ + - ]: 9291790 : auto alphamat = state[volfracIdx(nmat, imat)];
319 [ + - ]: 9291790 : auto arhomat = state[densityIdx(nmat, imat)];
320 : 9291790 : auto arhoemat = state[energyIdx(nmat, imat)];
321 : 9291790 : pri[pressureIdx(nmat,imat)] = eos_pressure< tag::multimat >(
322 [ + - ]: 9291790 : m_system, arhomat, vel[0], vel[1], vel[2], arhoemat, alphamat,
323 : : imat);
324 : 9291790 : pri[pressureIdx(nmat,imat)] = constrain_pressure< tag::multimat >(
325 : 9291790 : m_system, pri[pressureIdx(nmat,imat)], alphamat, imat);
326 : : }
327 : :
328 : : // Evaluate bulk velocity at quadrature point
329 [ + + ]: 16383208 : for (std::size_t idir=0; idir<3; ++idir) {
330 : 12287406 : pri[velocityIdx(nmat,idir)] = vel[idir];
331 : : }
332 : :
333 [ + + ]: 25674998 : for(std::size_t k = 0; k < nprim(); k++)
334 : : {
335 : 21579196 : auto mark = k * ndof;
336 [ + + ]: 21579196 : R[mark] += w * pri[k];
337 [ + + ]: 21579196 : if(ndof > 1)
338 : : {
339 [ + + ]: 46086400 : for(std::size_t idir = 0; idir < 3; idir++)
340 : 34564800 : R[mark+idir+1] += w * pri[k] * B[idir+1];
341 : : }
342 : : }
343 : : }
344 : :
345 : : // Update the DG solution of primitive variables
346 [ + + ]: 14614262 : for(std::size_t k = 0; k < nprim(); k++)
347 : : {
348 : 12361916 : auto mark = k * ndof;
349 : 12361916 : auto rmark = k * rdof;
350 [ + + ]: 12361916 : prim(e, rmark, 0) = R[mark] / L(e, mark, 0);
351 [ + + ]: 12361916 : if(ndof > 1)
352 : : {
353 [ + + ]: 9217280 : for(std::size_t idir = 0; idir < 3; idir++)
354 : : {
355 [ + + ]: 6912960 : prim(e, rmark+idir+1, 0) = R[mark+idir+1] / L(e, mark+idir+1, 0);
356 [ + + ]: 6912960 : if(fabs(prim(e, rmark+idir+1, 0)) < 1e-16)
357 : 2917061 : prim(e, rmark+idir+1, 0) = 0;
358 : : }
359 : : }
360 : : }
361 : : }
362 : 5315 : }
363 : :
364 : : //! Clean up the state of trace materials for this PDE system
365 : : //! \param[in] geoElem Element geometry array
366 : : //! \param[in,out] unk Array of unknowns
367 : : //! \param[in,out] prim Array of primitives
368 : : //! \param[in] nielem Number of internal elements
369 : : //! \details This function cleans up the state of materials present in trace
370 : : //! quantities in each cell. Specifically, the state of materials with
371 : : //! very low volume-fractions in a cell is replaced by the state of the
372 : : //! material which is present in the largest quantity in that cell. This
373 : : //! becomes necessary when shocks pass through cells which contain a very
374 : : //! small amount of material. The state of that tiny material might
375 : : //! become unphysical and cause solution to diverge; thus requiring such
376 : : //! a "reset".
377 : 5250 : void cleanTraceMaterial( const tk::Fields& geoElem,
378 : : tk::Fields& unk,
379 : : tk::Fields& prim,
380 : : std::size_t nielem ) const
381 : : {
382 : 5250 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
383 : 5250 : const auto nmat =
384 : 5250 : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[m_system];
385 : :
386 : : Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
387 : : "vector and primitive vector at recent time step incorrect" );
388 : : Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
389 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
390 : : Assert( prim.nprop() == rdof*nprim(), "Number of components in vector of "
391 : : "primitive quantities must equal "+ std::to_string(rdof*nprim()) );
392 : : Assert( (g_inputdeck.get< tag::discr, tag::ndof >()) <= 4, "High-order "
393 : : "discretizations not set up for multimat cleanTraceMaterial()" );
394 : :
395 : : auto al_eps = 1.0e-02;
396 : : auto neg_density = false;
397 : :
398 [ + + ]: 2235150 : for (std::size_t e=0; e<nielem; ++e)
399 : : {
400 : : // find material in largest quantity, and determine if pressure
401 : : // relaxation is needed. If it is, determine materials that need
402 : : // relaxation, and the total volume of these materials.
403 [ + - ]: 2229900 : std::vector< int > relaxInd(nmat, 0);
404 : : auto almax(0.0), relaxVol(0.0);
405 : : std::size_t kmax = 0;
406 [ + + ]: 7782600 : for (std::size_t k=0; k<nmat; ++k)
407 : : {
408 [ + + ]: 5552700 : auto al = unk(e, volfracDofIdx(nmat, k, rdof, 0), m_offset);
409 [ + + ]: 5552700 : if (al > almax)
410 : : {
411 : : almax = al;
412 : : kmax = k;
413 : : }
414 [ + + ]: 2112702 : else if (al < al_eps)
415 : : {
416 : 2101120 : relaxInd[k] = 1;
417 : : relaxVol += al;
418 : : }
419 : : }
420 [ + + ]: 2229900 : relaxInd[kmax] = 1;
421 : : relaxVol += almax;
422 : :
423 [ + + ]: 2229900 : auto u = prim(e, velocityDofIdx(nmat, 0, rdof, 0), m_offset);
424 : 2229900 : auto v = prim(e, velocityDofIdx(nmat, 1, rdof, 0), m_offset);
425 [ + + ]: 2229900 : auto w = prim(e, velocityDofIdx(nmat, 2, rdof, 0), m_offset);
426 [ + + ]: 2229900 : auto pmax = prim(e, pressureDofIdx(nmat, kmax, rdof, 0), m_offset)/almax;
427 [ + + ]: 2229900 : auto tmax = eos_temperature< eq >(m_system,
428 : : unk(e, densityDofIdx(nmat, kmax, rdof, 0), m_offset), u, v, w,
429 : : unk(e, energyDofIdx(nmat, kmax, rdof, 0), m_offset), almax, kmax);
430 : :
431 : : tk::real p_target(0.0), d_al(0.0), d_arE(0.0);
432 : : //// get equilibrium pressure
433 : : //std::vector< tk::real > kmat(nmat, 0.0);
434 : : //tk::real ratio(0.0);
435 : : //for (std::size_t k=0; k<nmat; ++k)
436 : : //{
437 : : // auto arhok = unk(e, densityDofIdx(nmat,k,rdof,0), m_offset);
438 : : // auto alk = unk(e, volfracDofIdx(nmat,k,rdof,0), m_offset);
439 : : // auto apk = prim(e, pressureDofIdx(nmat,k,rdof,0), m_offset);
440 : : // auto ak = eos_soundspeed< eq >(m_system, arhok, apk, alk, k );
441 : : // kmat[k] = arhok * ak * ak;
442 : :
443 : : // p_target += alk * apk / kmat[k];
444 : : // ratio += alk * alk / kmat[k];
445 : : //}
446 : : //p_target /= ratio;
447 : : //p_target = std::max(p_target, 1e-14);
448 [ + + ]: 2229931 : p_target = std::max(pmax, 1e-14);
449 : :
450 : : // 1. Correct minority materials and store volume/energy changes
451 [ + + ]: 7782600 : for (std::size_t k=0; k<nmat; ++k)
452 : : {
453 : 5552700 : auto alk = unk(e, volfracDofIdx(nmat, k, rdof, 0), m_offset);
454 [ + + ]: 5552700 : auto pk = prim(e, pressureDofIdx(nmat, k, rdof, 0), m_offset) / alk;
455 : : auto Pck = pstiff< eq >(m_system, k);
456 : : // for positive volume fractions
457 [ + + ]: 5552700 : if (matExists(alk))
458 : : {
459 : : // check if volume fraction is lesser than threshold (al_eps) and
460 : : // if the material (effective) pressure is negative. If either of
461 : : // these conditions is true, perform pressure relaxation.
462 [ + + ][ + + ]: 2422575 : if ((alk < al_eps) || (pk+Pck < 0.0)/*&& (std::fabs((pk-pmax)/pmax) > 1e-08)*/)
463 : : {
464 : : //auto gk = gamma< eq >(m_system, k);
465 : :
466 : : tk::real alk_new(0.0);
467 : : //// volume change based on polytropic expansion/isentropic compression
468 : : //if (pk > p_target)
469 : : //{
470 : : // alk_new = std::pow((pk/p_target), (1.0/gk)) * alk;
471 : : //}
472 : : //else
473 : : //{
474 : : // auto arhok = unk(e, densityDofIdx(nmat, k, rdof, 0), m_offset);
475 : : // auto ck = eos_soundspeed< eq >(m_system, arhok, alk*pk, alk, k);
476 : : // auto kk = arhok * ck * ck;
477 : : // alk_new = alk - (alk*alk/kk) * (p_target-pk);
478 : : //}
479 : : alk_new = alk;
480 : :
481 : : // energy change
482 : 157994 : auto rhomat = unk(e, densityDofIdx(nmat, k, rdof, 0), m_offset)
483 : 157994 : / alk_new;
484 : 157994 : auto rhoEmat = eos_totalenergy< eq >(m_system, rhomat, u, v, w,
485 : : p_target, k);
486 : :
487 : : // volume-fraction and total energy flux into majority material
488 : 157994 : d_al += (alk - alk_new);
489 : 157994 : d_arE += (unk(e, energyDofIdx(nmat, k, rdof, 0), m_offset)
490 : 157994 : - alk_new * rhoEmat);
491 : :
492 : : // update state of trace material
493 : : unk(e, volfracDofIdx(nmat, k, rdof, 0), m_offset) = alk_new;
494 : 157994 : unk(e, energyDofIdx(nmat, k, rdof, 0), m_offset) = alk_new*rhoEmat;
495 : 157994 : prim(e, pressureDofIdx(nmat, k, rdof, 0), m_offset) = alk_new*p_target;
496 : : }
497 : : }
498 : : // check for unbounded volume fractions
499 [ + + ]: 3130125 : else if (alk < 0.0)
500 : : {
501 : 14 : auto rhok = eos_density< eq >(m_system, p_target, tmax, k);
502 : 14 : d_al += (alk - 1e-14);
503 : : // update state of trace material
504 : 14 : unk(e, volfracDofIdx(nmat, k, rdof, 0), m_offset) = 1e-14;
505 : 14 : unk(e, densityDofIdx(nmat, k, rdof, 0), m_offset) = 1e-14 * rhok;
506 : 14 : unk(e, energyDofIdx(nmat, k, rdof, 0), m_offset) = 1e-14
507 : 14 : * eos_totalenergy< eq >(m_system, rhok, u, v, w, p_target, k);
508 : 14 : prim(e, pressureDofIdx(nmat, k, rdof, 0), m_offset) = 1e-14 *
509 : : p_target;
510 [ + + ]: 56 : for (std::size_t i=1; i<rdof; ++i) {
511 : 42 : unk(e, volfracDofIdx(nmat, k, rdof, i), m_offset) = 0.0;
512 : 42 : unk(e, densityDofIdx(nmat, k, rdof, i), m_offset) = 0.0;
513 : 42 : unk(e, energyDofIdx(nmat, k, rdof, i), m_offset) = 0.0;
514 : 42 : prim(e, pressureDofIdx(nmat, k, rdof, i), m_offset) = 0.0;
515 : : }
516 : : }
517 : : else {
518 : 3130111 : auto rhok = unk(e, densityDofIdx(nmat, k, rdof, 0), m_offset) / alk;
519 : : // update state of trace material
520 : 3130111 : unk(e, energyDofIdx(nmat, k, rdof, 0), m_offset) = alk
521 : 3130111 : * eos_totalenergy< eq >(m_system, rhok, u, v, w, p_target, k);
522 : 3130111 : prim(e, pressureDofIdx(nmat, k, rdof, 0), m_offset) = alk *
523 : : p_target;
524 [ + + ]: 5699878 : for (std::size_t i=1; i<rdof; ++i) {
525 : 2569767 : unk(e, energyDofIdx(nmat, k, rdof, i), m_offset) = 0.0;
526 : 2569767 : prim(e, pressureDofIdx(nmat, k, rdof, i), m_offset) = 0.0;
527 : : }
528 : : }
529 : : }
530 : :
531 : : // 2. Based on volume change in majority material, compute energy change
532 : : //auto gmax = gamma< eq >(m_system, kmax);
533 : : //auto pmax_new = pmax * std::pow(almax/(almax+d_al), gmax);
534 : : //auto rhomax_new = unk(e, densityDofIdx(nmat, kmax, rdof, 0), m_offset)
535 : : // / (almax+d_al);
536 : : //auto rhoEmax_new = eos_totalenergy< eq >(m_system, rhomax_new, u, v, w,
537 : : // pmax_new, kmax);
538 : : //auto d_arEmax_new = (almax+d_al) * rhoEmax_new
539 : : // - unk(e, energyDofIdx(nmat, kmax, rdof, 0), m_offset);
540 : :
541 : 2229900 : unk(e, volfracDofIdx(nmat, kmax, rdof, 0), m_offset) += d_al;
542 : : //unk(e, energyDofIdx(nmat, kmax, rdof, 0), m_offset) += d_arEmax_new;
543 : :
544 : : // 2. Flux energy change into majority material
545 : 2229900 : unk(e, energyDofIdx(nmat, kmax, rdof, 0), m_offset) += d_arE;
546 : 2229900 : prim(e, pressureDofIdx(nmat, kmax, rdof, 0), m_offset) =
547 [ + - ]: 2229900 : eos_pressure< eq >(m_system,
548 : : unk(e, densityDofIdx(nmat, kmax, rdof, 0), m_offset), u, v, w,
549 : : unk(e, energyDofIdx(nmat, kmax, rdof, 0), m_offset),
550 : : unk(e, volfracDofIdx(nmat, kmax, rdof, 0), m_offset), kmax);
551 : :
552 : : // enforce unit sum of volume fractions
553 : : auto alsum = 0.0;
554 [ + + ]: 7782600 : for (std::size_t k=0; k<nmat; ++k)
555 : 5552700 : alsum += unk(e, volfracDofIdx(nmat, k, rdof, 0), m_offset);
556 : :
557 [ + + ]: 7782600 : for (std::size_t k=0; k<nmat; ++k) {
558 : 5552700 : unk(e, volfracDofIdx(nmat, k, rdof, 0), m_offset) /= alsum;
559 : 5552700 : unk(e, densityDofIdx(nmat, k, rdof, 0), m_offset) /= alsum;
560 : 5552700 : unk(e, energyDofIdx(nmat, k, rdof, 0), m_offset) /= alsum;
561 : 5552700 : prim(e, pressureDofIdx(nmat, k, rdof, 0), m_offset) /= alsum;
562 : : }
563 : :
564 : : //// bulk quantities
565 : : //auto rhoEb(0.0), rhob(0.0), volb(0.0);
566 : : //for (std::size_t k=0; k<nmat; ++k)
567 : : //{
568 : : // if (relaxInd[k] > 0.0)
569 : : // {
570 : : // rhoEb += unk(e, energyDofIdx(nmat,k,rdof,0), m_offset);
571 : : // volb += unk(e, volfracDofIdx(nmat,k,rdof,0), m_offset);
572 : : // rhob += unk(e, densityDofIdx(nmat,k,rdof,0), m_offset);
573 : : // }
574 : : //}
575 : :
576 : : //// 2. find mixture-pressure
577 : : //tk::real pmix(0.0), den(0.0);
578 : : //pmix = rhoEb - 0.5*rhob*(u*u+v*v+w*w);
579 : : //for (std::size_t k=0; k<nmat; ++k)
580 : : //{
581 : : // auto gk = gamma< eq >(m_system, k);
582 : : // auto Pck = pstiff< eq >(m_system, k);
583 : :
584 : : // pmix -= unk(e, volfracDofIdx(nmat,k,rdof,0), m_offset) * gk * Pck *
585 : : // relaxInd[k] / (gk-1.0);
586 : : // den += unk(e, volfracDofIdx(nmat,k,rdof,0), m_offset) * relaxInd[k]
587 : : // / (gk-1.0);
588 : : //}
589 : : //pmix /= den;
590 : :
591 : : //// 3. correct energies
592 : : //for (std::size_t k=0; k<nmat; ++k)
593 : : //{
594 : : // if (relaxInd[k] > 0.0)
595 : : // {
596 : : // auto alk_new = unk(e, volfracDofIdx(nmat,k,rdof,0), m_offset);
597 : : // unk(e, energyDofIdx(nmat,k,rdof,0), m_offset) = alk_new *
598 : : // eos_totalenergy< eq >(m_system, rhomat[k], u, v, w, pmix, k);
599 : : // prim(e, pressureDofIdx(nmat, k, rdof, 0), m_offset) = alk_new * pmix;
600 : : // }
601 : : //}
602 : :
603 : 2229900 : pmax = prim(e, pressureDofIdx(nmat, kmax, rdof, 0), m_offset) /
604 : : unk(e, volfracDofIdx(nmat, kmax, rdof, 0), m_offset);
605 : :
606 : : // check for unphysical state
607 [ + + ]: 7782600 : for (std::size_t k=0; k<nmat; ++k)
608 : : {
609 [ - + ]: 5552700 : auto alpha = unk(e, volfracDofIdx(nmat, k, rdof, 0), m_offset);
610 : 5552700 : auto arho = unk(e, densityDofIdx(nmat, k, rdof, 0), m_offset);
611 : 5552700 : auto apr = prim(e, pressureDofIdx(nmat, k, rdof, 0), m_offset);
612 : :
613 : : // lambda for screen outputs
614 : 5552700 : auto screenout = [&]()
615 : : {
616 : 0 : std::cout << "Element centroid: " << geoElem(e,1,0) << ", "
617 : 0 : << geoElem(e,2,0) << ", " << geoElem(e,3,0) << std::endl;
618 : 0 : std::cout << "Material-id: " << k << std::endl;
619 : 0 : std::cout << "Volume-fraction: " << alpha << std::endl;
620 : 0 : std::cout << "Partial density: " << arho << std::endl;
621 : 0 : std::cout << "Partial pressure: " << apr << std::endl;
622 : 0 : std::cout << "Major pressure: " << pmax << std::endl;
623 : 0 : std::cout << "Major temperature:" << tmax << std::endl;
624 : 0 : std::cout << "Velocity: " << u << ", " << v << ", " << w
625 : : << std::endl;
626 : : };
627 : :
628 [ - + ]: 5552700 : if (arho < 0.0)
629 : : {
630 : : neg_density = true;
631 [ - - ]: 0 : screenout();
632 : : }
633 : : }
634 : : }
635 : :
636 [ - + ][ - - ]: 5250 : if (neg_density) Throw("Negative partial density.");
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
637 : 5250 : }
638 : :
639 : : //! Reconstruct second-order solution from first-order
640 : : //! \param[in] geoElem Element geometry array
641 : : //! \param[in] fd Face connectivity and boundary conditions object
642 : : //! \param[in] esup Elements-surrounding-nodes connectivity
643 : : //! \param[in] inpoel Element-node connectivity
644 : : //! \param[in] coord Array of nodal coordinates
645 : : //! \param[in,out] U Solution vector at recent time step
646 : : //! \param[in,out] P Vector of primitives at recent time step
647 : 4125 : void reconstruct( tk::real,
648 : : const tk::Fields&,
649 : : const tk::Fields& geoElem,
650 : : const inciter::FaceData& fd,
651 : : const std::map< std::size_t, std::vector< std::size_t > >&
652 : : esup,
653 : : const std::vector< std::size_t >& inpoel,
654 : : const tk::UnsMesh::Coords& coord,
655 : : tk::Fields& U,
656 : : tk::Fields& P ) const
657 : : {
658 : 4125 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
659 : 4125 : const auto intsharp =
660 [ + - ]: 4125 : g_inputdeck.get< tag::param, tag::multimat, tag::intsharp >()[m_system];
661 : :
662 : : bool is_p0p1(false);
663 [ + - ][ + + ]: 4125 : if (rdof == 4 && g_inputdeck.get< tag::discr, tag::ndof >() == 1)
664 : : is_p0p1 = true;
665 : :
666 : : // do reconstruction only if P0P1 or if interface reconstruction is active
667 [ + + ]: 4125 : if (is_p0p1 || (intsharp > 0)) {
668 [ + + ]: 3750 : const auto nelem = fd.Esuel().size()/4;
669 [ + + ]: 3750 : const auto nmat =
670 : : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[m_system];
671 : :
672 : : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
673 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
674 : :
675 : : //----- reconstruction of conserved quantities -----
676 : : //--------------------------------------------------
677 : : // specify how many variables need to be reconstructed
678 : 3750 : std::array< std::size_t, 2 > varRange {{0, m_ncomp-1}};
679 [ + + ]: 3750 : if (!is_p0p1)
680 : 375 : varRange = {{volfracIdx(nmat, 0), volfracIdx(nmat, nmat-1)}};
681 : :
682 : : // 1. solve 3x3 least-squares system
683 [ + + ]: 685950 : for (std::size_t e=0; e<nelem; ++e)
684 : : {
685 : : // Reconstruct second-order dofs of volume-fractions in Taylor space
686 : : // using nodal-stencils, for a good interface-normal estimate
687 : 682200 : tk::recoLeastSqExtStencil( rdof, m_offset, e, esup, inpoel, geoElem,
688 : : U, varRange );
689 : : }
690 : :
691 : : // 2. transform reconstructed derivatives to Dubiner dofs
692 : 3750 : tk::transform_P0P1(m_offset, rdof, nelem, inpoel, coord, U, varRange);
693 : :
694 : : //----- reconstruction of primitive quantities -----
695 : : //--------------------------------------------------
696 : : // For multimat, conserved and primitive quantities are reconstructed
697 : : // separately.
698 [ + + ]: 3750 : if (is_p0p1) {
699 : : // 1.
700 [ + + ]: 458175 : for (std::size_t e=0; e<nelem; ++e)
701 : : {
702 : : // Reconstruct second-order dofs of volume-fractions in Taylor space
703 : : // using nodal-stencils, for a good interface-normal estimate
704 : 454800 : tk::recoLeastSqExtStencil( rdof, m_offset, e, esup, inpoel, geoElem,
705 : 454800 : P, {0, nprim()-1} );
706 : : }
707 : :
708 : : // 2.
709 : 3375 : tk::transform_P0P1(m_offset, rdof, nelem, inpoel, coord, P,
710 : 3375 : {0, nprim()-1});
711 : : }
712 : : }
713 : 4125 : }
714 : :
715 : : //! Limit second-order solution, and primitive quantities separately
716 : : //! \param[in] t Physical time
717 : : //! \param[in] geoFace Face geometry array
718 : : //! \param[in] geoElem Element geometry array
719 : : //! \param[in] fd Face connectivity and boundary conditions object
720 : : //! \param[in] esup Elements-surrounding-nodes connectivity
721 : : //! \param[in] inpoel Element-node connectivity
722 : : //! \param[in] coord Array of nodal coordinates
723 : : //! \param[in] ndofel Vector of local number of degrees of freedome
724 : : // //! \param[in] gid Local->global node id map
725 : : // //! \param[in] bid Local chare-boundary node ids (value) associated to
726 : : // //! global node ids (key)
727 : : // //! \param[in] uNodalExtrm Chare-boundary nodal extrema for conservative
728 : : // //! variables
729 : : // //! \param[in] pNodalExtrm Chare-boundary nodal extrema for primitive
730 : : // //! variables
731 : : //! \param[in,out] U Solution vector at recent time step
732 : : //! \param[in,out] P Vector of primitives at recent time step
733 : 4125 : void limit( [[maybe_unused]] tk::real t,
734 : : const tk::Fields& geoFace,
735 : : const tk::Fields& geoElem,
736 : : const inciter::FaceData& fd,
737 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
738 : : const std::vector< std::size_t >& inpoel,
739 : : const tk::UnsMesh::Coords& coord,
740 : : const std::vector< std::size_t >& ndofel,
741 : : const std::vector< std::size_t >&,
742 : : const std::unordered_map< std::size_t, std::size_t >&,
743 : : const std::vector< std::vector<tk::real> >&,
744 : : const std::vector< std::vector<tk::real> >&,
745 : : tk::Fields& U,
746 : : tk::Fields& P,
747 : : std::vector< std::size_t >& shockmarker ) const
748 : : {
749 : : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
750 : : "vector and primitive vector at recent time step incorrect" );
751 : :
752 : 4125 : const auto limiter = g_inputdeck.get< tag::discr, tag::limiter >();
753 : 4125 : const auto nmat =
754 [ - + ]: 4125 : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[m_system];
755 : :
756 : : // limit vectors of conserved and primitive quantities
757 [ - + ]: 4125 : if (limiter == ctr::LimiterType::SUPERBEEP1)
758 : : {
759 : 0 : SuperbeeMultiMat_P1( fd.Esuel(), inpoel, ndofel, m_system, m_offset,
760 : : coord, U, P, nmat );
761 : : }
762 [ + - ]: 4125 : else if (limiter == ctr::LimiterType::VERTEXBASEDP1)
763 : : {
764 : 4125 : VertexBasedMultiMat_P1( esup, inpoel, ndofel, fd.Esuel().size()/4,
765 : 4125 : m_system, m_offset, fd, geoFace, geoElem, coord, U, P, nmat,
766 : : shockmarker );
767 : : }
768 : : else
769 : : {
770 [ - - ][ - - ]: 0 : Throw("Limiter type not configured for multimat.");
[ - - ][ - - ]
[ - - ][ - - ]
771 : : }
772 : 4125 : }
773 : :
774 : : //! Compute right hand side
775 : : //! \param[in] t Physical time
776 : : //! \param[in] geoFace Face geometry array
777 : : //! \param[in] geoElem Element geometry array
778 : : //! \param[in] fd Face connectivity and boundary conditions object
779 : : //! \param[in] inpoel Element-node connectivity
780 : : //! \param[in] coord Array of nodal coordinates
781 : : //! \param[in] U Solution vector at recent time step
782 : : //! \param[in] P Primitive vector at recent time step
783 : : //! \param[in] ndofel Vector of local number of degrees of freedome
784 : : //! \param[in,out] R Right-hand side vector computed
785 : 5250 : void rhs( tk::real t,
786 : : const tk::Fields& geoFace,
787 : : const tk::Fields& geoElem,
788 : : const inciter::FaceData& fd,
789 : : const std::vector< std::size_t >& inpoel,
790 : : const std::vector< std::unordered_set< std::size_t > >&,
791 : : const tk::UnsMesh::Coords& coord,
792 : : const tk::Fields& U,
793 : : const tk::Fields& P,
794 : : const std::vector< std::size_t >& ndofel,
795 : : tk::Fields& R ) const
796 : : {
797 : 5250 : const auto ndof = g_inputdeck.get< tag::discr, tag::ndof >();
798 : 5250 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
799 : 5250 : const auto nmat =
800 : 5250 : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[m_system];
801 : 5250 : const auto intsharp =
802 : : g_inputdeck.get< tag::param, tag::multimat, tag::intsharp >()[m_system];
803 : :
804 : 5250 : const auto nelem = fd.Esuel().size()/4;
805 : :
806 : : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
807 : : "vector and primitive vector at recent time step incorrect" );
808 : : Assert( U.nunk() == R.nunk(), "Number of unknowns in solution "
809 : : "vector and right-hand side at recent time step incorrect" );
810 : : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
811 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
812 : : Assert( P.nprop() == rdof*nprim(), "Number of components in primitive "
813 : : "vector must equal "+ std::to_string(rdof*nprim()) );
814 : : Assert( R.nprop() == ndof*m_ncomp, "Number of components in right-hand "
815 : : "side vector must equal "+ std::to_string(ndof*m_ncomp) );
816 : : Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
817 : : "Mismatch in inpofa size" );
818 : : Assert( ndof <= 4, "DGP2 not set up for multi-material" );
819 : :
820 : : // set rhs to zero
821 : : R.fill(0.0);
822 : :
823 : : // allocate space for Riemann derivatives used in non-conservative terms
824 : : std::vector< std::vector< tk::real > >
825 [ + - ][ + - ]: 10500 : riemannDeriv( 3*nmat+1, std::vector<tk::real>(U.nunk(),0.0) );
826 : :
827 : : // vectors to store the data of riemann velocity used for reconstruction
828 : : // in volume fraction equation
829 [ + - ]: 10500 : std::vector< std::vector< tk::real > > vriem( U.nunk() );
830 [ + - ][ + - ]: 10500 : std::vector< std::vector< tk::real > > riemannLoc( U.nunk() );
831 : :
832 : : // configure Riemann flux function
833 : : auto rieflxfn =
834 : : [this]( const std::array< tk::real, 3 >& fn,
835 : : const std::array< std::vector< tk::real >, 2 >& u,
836 : : const std::vector< std::array< tk::real, 3 > >& v )
837 : : { return m_riemann.flux( fn, u, v ); };
838 : :
839 : : // configure a no-op lambda for prescribed velocity
840 : : auto velfn = [this]( ncomp_t, ncomp_t, tk::real, tk::real, tk::real,
841 : : tk::real ){
842 : 9640575 : return std::vector< std::array< tk::real, 3 > >( m_ncomp ); };
843 : :
844 : : // compute internal surface flux integrals
845 [ + - ][ + - ]: 10500 : tk::surfInt( m_system, nmat, m_offset, t, ndof, rdof, inpoel, coord,
[ - - ]
846 : : fd, geoFace, geoElem, rieflxfn, velfn, U, P, ndofel, R,
847 : : vriem, riemannLoc, riemannDeriv, intsharp );
848 : :
849 [ + + ]: 5250 : if(ndof > 1)
850 : : // compute volume integrals
851 [ + - ][ + - ]: 2250 : tk::volInt( m_system, nmat, m_offset, t, ndof, rdof, nelem, inpoel,
[ - - ]
852 : : coord, geoElem, flux, velfn, U, P, ndofel, R, intsharp );
853 : :
854 : : // compute boundary surface flux integrals
855 [ + + ]: 36750 : for (const auto& b : m_bc)
856 [ + - ][ + - ]: 94500 : tk::bndSurfInt( m_system, nmat, m_offset, ndof, rdof, b.first,
[ - - ]
857 : : fd, geoFace, geoElem, inpoel, coord, t, rieflxfn, velfn,
858 [ + - ]: 31500 : b.second, U, P, ndofel, R, vriem, riemannLoc,
859 : : riemannDeriv, intsharp );
860 : :
861 : : Assert( riemannDeriv.size() == 3*nmat+1, "Size of Riemann derivative "
862 : : "vector incorrect" );
863 : :
864 : : // get derivatives from riemannDeriv
865 [ + + ]: 44250 : for (std::size_t k=0; k<riemannDeriv.size(); ++k)
866 : : {
867 : : Assert( riemannDeriv[k].size() == U.nunk(), "Riemann derivative vector "
868 : : "for non-conservative terms has incorrect size" );
869 [ + + ]: 26730000 : for (std::size_t e=0; e<U.nunk(); ++e)
870 : 26691000 : riemannDeriv[k][e] /= geoElem(e, 0, 0);
871 : : }
872 : :
873 : : std::vector< std::vector< tk::real > >
874 [ + - ][ + - ]: 10500 : vriempoly( U.nunk(), std::vector<tk::real>(12,0.0) );
875 : : // get the polynomial solution of Riemann velocity at the interface.
876 : : // not required if interface reconstruction is used, since then volfrac
877 : : // equation is discretized using p0p1.
878 [ + + ]: 5250 : if (ndof > 1 && intsharp == 0)
879 [ + - ]: 750 : vriempoly = tk::solvevriem(nelem, vriem, riemannLoc);
880 : :
881 : : // compute volume integrals of non-conservative terms
882 [ + - ]: 5250 : tk::nonConservativeInt( m_system, nmat, m_offset, ndof, rdof, nelem,
883 : : inpoel, coord, geoElem, U, P, riemannDeriv,
884 : : vriempoly, ndofel, R, intsharp );
885 : :
886 : : // compute finite pressure relaxation terms
887 [ + + ]: 5250 : if (g_inputdeck.get< tag::param, tag::multimat, tag::prelax >()[m_system])
888 : : {
889 [ + - ]: 450 : const auto ct = g_inputdeck.get< tag::param, tag::multimat,
890 : : tag::prelax_timescale >()[m_system];
891 [ + - ]: 450 : tk::pressureRelaxationInt( m_system, nmat, m_offset, ndof, rdof, nelem,
892 : : inpoel, coord, geoElem, U, P, ndofel, ct, R,
893 : : intsharp );
894 : : }
895 : 5250 : }
896 : :
897 : : //! Evaluate the adaptive indicator and mark the ndof for each element
898 : : //! \param[in] nunk Number of unknowns
899 : : //! \param[in] coord Array of nodal coordinates
900 : : //! \param[in] inpoel Element-node connectivity
901 : : //! \param[in] fd Face connectivity and boundary conditions object
902 : : //! \param[in] unk Array of unknowns
903 : : //! \param[in] indicator p-refinement indicator type
904 : : //! \param[in] ndof Number of degrees of freedom in the solution
905 : : //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
906 : : //! \param[in] tolref Tolerance for p-refinement
907 : : //! \param[in,out] ndofel Vector of local number of degrees of freedome
908 [ - - ]: 0 : void eval_ndof( std::size_t nunk,
909 : : [[maybe_unused]] const tk::UnsMesh::Coords& coord,
910 : : [[maybe_unused]] const std::vector< std::size_t >& inpoel,
911 : : const inciter::FaceData& fd,
912 : : const tk::Fields& unk,
913 : : inciter::ctr::PrefIndicatorType indicator,
914 : : std::size_t ndof,
915 : : std::size_t ndofmax,
916 : : tk::real tolref,
917 : : std::vector< std::size_t >& ndofel ) const
918 : : {
919 : : const auto& esuel = fd.Esuel();
920 : 0 : const auto nmat =
921 [ - - ]: 0 : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[m_system];
922 : :
923 [ - - ]: 0 : if(indicator == inciter::ctr::PrefIndicatorType::SPECTRAL_DECAY)
924 : 0 : spectral_decay(nmat, nunk, esuel, unk, ndof, ndofmax, tolref, ndofel);
925 : : else
926 [ - - ][ - - ]: 0 : Throw( "No such adaptive indicator type" );
[ - - ][ - - ]
[ - - ][ - - ]
927 : 0 : }
928 : :
929 : : //! Compute the minimum time step size
930 : : //! \param[in] fd Face connectivity and boundary conditions object
931 : : //! \param[in] geoFace Face geometry array
932 : : //! \param[in] geoElem Element geometry array
933 : : // //! \param[in] ndofel Vector of local number of degrees of freedom
934 : : //! \param[in] U Solution vector at recent time step
935 : : //! \param[in] P Vector of primitive quantities at recent time step
936 : : //! \param[in] nielem Number of internal elements
937 : : //! \return Minimum time step size
938 : : //! \details The allowable dt is calculated by looking at the maximum
939 : : //! wave-speed in elements surrounding each face, times the area of that
940 : : //! face. Once the maximum of this quantity over the mesh is determined,
941 : : //! the volume of each cell is divided by this quantity. A minimum of this
942 : : //! ratio is found over the entire mesh, which gives the allowable dt.
943 : 275 : tk::real dt( const std::array< std::vector< tk::real >, 3 >&,
944 : : const std::vector< std::size_t >&,
945 : : const inciter::FaceData& fd,
946 : : const tk::Fields& geoFace,
947 : : const tk::Fields& geoElem,
948 : : const std::vector< std::size_t >& /*ndofel*/,
949 : : const tk::Fields& U,
950 : : const tk::Fields& P,
951 : : const std::size_t nielem ) const
952 : : {
953 : 275 : const auto ndof = g_inputdeck.get< tag::discr, tag::ndof >();
954 : 275 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
955 : 275 : const auto nmat =
956 : 275 : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[m_system];
957 : :
958 : : const auto& esuf = fd.Esuf();
959 : :
960 : : tk::real u, v, w, a, vn, dSV_l, dSV_r;
961 : 275 : std::vector< tk::real > delt(U.nunk(), 0.0);
962 [ + - ][ + - ]: 275 : std::vector< tk::real > ugp(m_ncomp, 0.0), pgp(nprim(), 0.0);
[ - - ][ - - ]
963 : :
964 : : // compute maximum characteristic speed at all internal element faces
965 [ + + ]: 436325 : for (std::size_t f=0; f<esuf.size()/2; ++f)
966 : : {
967 [ + - ]: 436050 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
968 : 436050 : auto er = esuf[2*f+1];
969 : :
970 : : // left element
971 : :
972 : : // Compute the basis function for the left element
973 [ + - ][ - - ]: 436050 : std::vector< tk::real > B_l(rdof, 0.0);
974 : 436050 : B_l[0] = 1.0;
975 : :
976 : : // get conserved quantities
977 [ + - ]: 436050 : ugp = eval_state( m_ncomp, m_offset, rdof, ndof, el, U, B_l, {0, m_ncomp-1} );
978 : : // get primitive quantities
979 [ + - ]: 436050 : pgp = eval_state( nprim(), m_offset, rdof, ndof, el, P, B_l, {0, nprim()-1} );
980 : :
981 : : // advection velocity
982 : 436050 : u = pgp[velocityIdx(nmat, 0)];
983 : 436050 : v = pgp[velocityIdx(nmat, 1)];
984 : 436050 : w = pgp[velocityIdx(nmat, 2)];
985 : :
986 : 436050 : vn = u*geoFace(f,1,0) + v*geoFace(f,2,0) + w*geoFace(f,3,0);
987 : :
988 : : // acoustic speed
989 : 436050 : a = 0.0;
990 [ + + ]: 1308150 : for (std::size_t k=0; k<nmat; ++k)
991 : : {
992 [ + + ]: 872100 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04) {
993 [ + - ][ + + ]: 895997 : a = std::max( a, eos_soundspeed< tag::multimat >( 0,
994 : : ugp[densityIdx(nmat, k)], pgp[pressureIdx(nmat, k)],
995 : : ugp[volfracIdx(nmat, k)], k ) );
996 : : }
997 : : }
998 : :
999 : 436050 : dSV_l = geoFace(f,0,0) * (std::fabs(vn) + a);
1000 : :
1001 : : // right element
1002 [ + + ]: 436050 : if (er > -1) {
1003 : 336550 : std::size_t eR = static_cast< std::size_t >( er );
1004 : :
1005 : : // Compute the basis function for the right element
1006 [ + - ][ - - ]: 336550 : std::vector< tk::real > B_r(rdof, 0.0);
1007 : 336550 : B_r[0] = 1.0;
1008 : :
1009 : : // get conserved quantities
1010 [ + - ]: 336550 : ugp = eval_state( m_ncomp, m_offset, rdof, ndof, eR, U, B_r, {0, m_ncomp-1});
1011 : : // get primitive quantities
1012 [ + - ][ - - ]: 336550 : pgp = eval_state( nprim(), m_offset, rdof, ndof, eR, P, B_r, {0, nprim()-1});
1013 : :
1014 : : // advection velocity
1015 : 336550 : u = pgp[velocityIdx(nmat, 0)];
1016 : 336550 : v = pgp[velocityIdx(nmat, 1)];
1017 : 336550 : w = pgp[velocityIdx(nmat, 2)];
1018 : :
1019 : 336550 : vn = u*geoFace(f,1,0) + v*geoFace(f,2,0) + w*geoFace(f,3,0);
1020 : :
1021 : : // acoustic speed
1022 : 336550 : a = 0.0;
1023 [ + + ]: 1009650 : for (std::size_t k=0; k<nmat; ++k)
1024 : : {
1025 [ + + ]: 673100 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04) {
1026 [ + - ][ + + ]: 691762 : a = std::max( a, eos_soundspeed< tag::multimat >( 0,
1027 : : ugp[densityIdx(nmat, k)], pgp[pressureIdx(nmat, k)],
1028 : : ugp[volfracIdx(nmat, k)], k ) );
1029 : : }
1030 : : }
1031 : :
1032 [ + + ]: 336550 : dSV_r = geoFace(f,0,0) * (std::fabs(vn) + a);
1033 : :
1034 [ + - ]: 336550 : delt[eR] += std::max( dSV_l, dSV_r );
1035 : : } else {
1036 : 99500 : dSV_r = dSV_l;
1037 : : }
1038 : :
1039 [ + - ]: 436050 : delt[el] += std::max( dSV_l, dSV_r );
1040 : : }
1041 : :
1042 : 275 : tk::real mindt = std::numeric_limits< tk::real >::max();
1043 : :
1044 : : // compute allowable dt
1045 [ + + ]: 189775 : for (std::size_t e=0; e<nielem; ++e)
1046 : : {
1047 [ + + ]: 191383 : mindt = std::min( mindt, geoElem(e,0,0)/delt[e] );
1048 : : }
1049 : :
1050 : : tk::real dgp = 0.0;
1051 [ + + ]: 275 : if (ndof == 4)
1052 : : {
1053 : : dgp = 1.0;
1054 : : }
1055 [ - + ]: 25 : else if (ndof == 10)
1056 : : {
1057 : : dgp = 2.0;
1058 : : }
1059 : :
1060 : : // Scale smallest dt with CFL coefficient and the CFL is scaled by (2*p+1)
1061 : : // where p is the order of the DG polynomial by linear stability theory.
1062 [ + - ]: 275 : mindt /= (2.0*dgp + 1.0);
1063 : 275 : return mindt;
1064 : : }
1065 : :
1066 : : //! Extract the velocity field at cell nodes. Currently unused.
1067 : : //! \param[in] U Solution vector at recent time step
1068 : : //! \param[in] N Element node indices
1069 : : //! \return Array of the four values of the velocity field
1070 : : std::array< std::array< tk::real, 4 >, 3 >
1071 : : velocity( const tk::Fields& U,
1072 : : const std::array< std::vector< tk::real >, 3 >&,
1073 : : const std::array< std::size_t, 4 >& N ) const
1074 : : {
1075 : : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
1076 : : const auto nmat =
1077 : : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[0];
1078 : :
1079 : : std::array< std::array< tk::real, 4 >, 3 > v;
1080 : : v[0] = U.extract( momentumDofIdx(nmat, 0, rdof, 0), m_offset, N );
1081 : : v[1] = U.extract( momentumDofIdx(nmat, 1, rdof, 0), m_offset, N );
1082 : : v[2] = U.extract( momentumDofIdx(nmat, 2, rdof, 0), m_offset, N );
1083 : :
1084 : : std::vector< std::array< tk::real, 4 > > ar;
1085 : : ar.resize(nmat);
1086 : : for (std::size_t k=0; k<nmat; ++k)
1087 : : ar[k] = U.extract( densityDofIdx(nmat, k, rdof, 0), m_offset, N );
1088 : :
1089 : : std::array< tk::real, 4 > r{{ 0.0, 0.0, 0.0, 0.0 }};
1090 : : for (std::size_t i=0; i<r.size(); ++i) {
1091 : : for (std::size_t k=0; k<nmat; ++k)
1092 : : r[i] += ar[k][i];
1093 : : }
1094 : :
1095 : : std::transform( r.begin(), r.end(), v[0].begin(), v[0].begin(),
1096 : : []( tk::real s, tk::real& d ){ return d /= s; } );
1097 : : std::transform( r.begin(), r.end(), v[1].begin(), v[1].begin(),
1098 : : []( tk::real s, tk::real& d ){ return d /= s; } );
1099 : : std::transform( r.begin(), r.end(), v[2].begin(), v[2].begin(),
1100 : : []( tk::real s, tk::real& d ){ return d /= s; } );
1101 : : return v;
1102 : : }
1103 : :
1104 : : //! Return analytic field names to be output to file
1105 : : //! \return Vector of strings labelling analytic fields output in file
1106 : : std::vector< std::string > analyticFieldNames() const {
1107 : 0 : auto nmat =
1108 : : g_inputdeck.get< tag::param, eq, tag::nmat >()[m_system];
1109 : :
1110 : 0 : return MultiMatFieldNames(nmat);
1111 : : }
1112 : :
1113 : : //! Return field names to be output to file
1114 : : //! \return Vector of strings labelling fields output in file
1115 : : std::vector< std::string > nodalFieldNames() const
1116 : : {
1117 : : auto nmat =
1118 : : g_inputdeck.get< tag::param, eq, tag::nmat >()[m_system];
1119 : :
1120 : : return MultiMatFieldNames(nmat);
1121 : : }
1122 : :
1123 : : //! Return time history field names to be output to file
1124 : : //! \return Vector of strings labelling time history fields output in file
1125 : : std::vector< std::string > histNames() const {
1126 : 0 : return MultiMatHistNames();
1127 : : }
1128 : :
1129 : : //! Return surface field output going to file
1130 : : std::vector< std::vector< tk::real > >
1131 : : surfOutput( const std::map< int, std::vector< std::size_t > >&,
1132 : : tk::Fields& ) const
1133 : : {
1134 : : std::vector< std::vector< tk::real > > s; // punt for now
1135 : : return s;
1136 : : }
1137 : :
1138 : : //! Return time history field output evaluated at time history points
1139 : : //! \param[in] h History point data
1140 : : //! \param[in] inpoel Element-node connectivity
1141 : : //! \param[in] coord Array of nodal coordinates
1142 : : //! \param[in] U Array of unknowns
1143 : : //! \param[in] P Array of primitive quantities
1144 : : //! \return Vector of time history output of bulk flow quantities (density,
1145 : : //! velocity, total energy, and pressure) evaluated at time history points
1146 : : std::vector< std::vector< tk::real > >
1147 : 0 : histOutput( const std::vector< HistData >& h,
1148 : : const std::vector< std::size_t >& inpoel,
1149 : : const tk::UnsMesh::Coords& coord,
1150 : : const tk::Fields& U,
1151 : : const tk::Fields& P ) const
1152 : : {
1153 : 0 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
1154 : 0 : const auto nmat =
1155 : 0 : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[m_system];
1156 : :
1157 : : const auto& x = coord[0];
1158 : : const auto& y = coord[1];
1159 : : const auto& z = coord[2];
1160 : :
1161 : 0 : std::vector< std::vector< tk::real > > Up(h.size());
1162 : :
1163 : : std::size_t j = 0;
1164 [ - - ]: 0 : for (const auto& p : h) {
1165 : 0 : auto e = p.get< tag::elem >();
1166 : 0 : auto chp = p.get< tag::coord >();
1167 : :
1168 : : // Evaluate inverse Jacobian
1169 : 0 : std::array< std::array< tk::real, 3>, 4 > cp{{
1170 : : {{ x[inpoel[4*e ]], y[inpoel[4*e ]], z[inpoel[4*e ]] }},
1171 : : {{ x[inpoel[4*e+1]], y[inpoel[4*e+1]], z[inpoel[4*e+1]] }},
1172 : : {{ x[inpoel[4*e+2]], y[inpoel[4*e+2]], z[inpoel[4*e+2]] }},
1173 : : {{ x[inpoel[4*e+3]], y[inpoel[4*e+3]], z[inpoel[4*e+3]] }} }};
1174 : 0 : auto J = tk::inverseJacobian( cp[0], cp[1], cp[2], cp[3] );
1175 : :
1176 : : // evaluate solution at history-point
1177 : 0 : std::array< tk::real, 3 > dc{{chp[0]-cp[0][0], chp[1]-cp[0][1],
1178 [ - - ]: 0 : chp[2]-cp[0][2]}};
1179 [ - - ]: 0 : auto B = tk::eval_basis(rdof, tk::dot(J[0],dc), tk::dot(J[1],dc),
1180 : : tk::dot(J[2],dc));
1181 [ - - ][ - - ]: 0 : auto uhp = eval_state(m_ncomp, m_offset, rdof, rdof, e, U, B, {0, m_ncomp-1});
1182 [ - - ][ - - ]: 0 : auto php = eval_state(nprim(), m_offset, rdof, rdof, e, P, B, {0, nprim()-1});
1183 : :
1184 : : // store solution in history output vector
1185 [ - - ][ - - ]: 0 : Up[j].resize(6, 0.0);
1186 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
1187 : 0 : Up[j][0] += uhp[densityIdx(nmat,k)];
1188 : 0 : Up[j][4] += uhp[energyIdx(nmat,k)];
1189 : 0 : Up[j][5] += php[pressureIdx(nmat,k)];
1190 : : }
1191 [ - - ]: 0 : Up[j][1] = php[velocityIdx(nmat,0)];
1192 : 0 : Up[j][2] = php[velocityIdx(nmat,1)];
1193 : 0 : Up[j][3] = php[velocityIdx(nmat,2)];
1194 [ - - ]: 0 : ++j;
1195 : : }
1196 : :
1197 : 0 : return Up;
1198 : : }
1199 : :
1200 : : //! Return names of integral variables to be output to diagnostics file
1201 : : //! \return Vector of strings labelling integral variables output
1202 : : std::vector< std::string > names() const
1203 : 12 : { return Problem::names( m_ncomp ); }
1204 : :
1205 : : //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
1206 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
1207 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
1208 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
1209 : : //! \param[in] t Physical time at which to evaluate the analytic solution
1210 : : //! \return Vector of analytic solution at given location and time
1211 : : std::vector< tk::real >
1212 : : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
1213 : : { return Problem::analyticSolution( m_system, m_ncomp, xi, yi, zi, t ); }
1214 : :
1215 : : //! Return analytic solution for conserved variables
1216 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
1217 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
1218 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
1219 : : //! \param[in] t Physical time at which to evaluate the analytic solution
1220 : : //! \return Vector of analytic solution at given location and time
1221 : : std::vector< tk::real >
1222 : : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
1223 : 1198100 : { return Problem::initialize( m_system, m_ncomp, xi, yi, zi, t ); }
1224 : :
1225 : : private:
1226 : : //! Equation system index
1227 : : const ncomp_t m_system;
1228 : : //! Number of components in this PDE system
1229 : : const ncomp_t m_ncomp;
1230 : : //! Offset PDE system operates from
1231 : : const ncomp_t m_offset;
1232 : : //! Riemann solver
1233 : : RiemannSolver m_riemann;
1234 : : //! BC configuration
1235 : : BCStateFn m_bc;
1236 : :
1237 : : //! Evaluate conservative part of physical flux function for this PDE system
1238 : : //! \param[in] system Equation system index
1239 : : //! \param[in] ncomp Number of scalar components in this PDE system
1240 : : //! \param[in] ugp Numerical solution at the Gauss point at which to
1241 : : //! evaluate the flux
1242 : : //! \return Flux vectors for all components in this PDE system
1243 : : //! \note The function signature must follow tk::FluxFn
1244 : : static tk::FluxFn::result_type
1245 : 2274000 : flux( ncomp_t system,
1246 : : [[maybe_unused]] ncomp_t ncomp,
1247 : : const std::vector< tk::real >& ugp,
1248 : : const std::vector< std::array< tk::real, 3 > >& )
1249 : : {
1250 : : //Assert( ugp.size() == ncomp, "Size mismatch" );
1251 : 2274000 : const auto nmat =
1252 : : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[system];
1253 : :
1254 : : tk::real rho(0.0), p(0.0);
1255 : : for (std::size_t k=0; k<nmat; ++k)
1256 : : rho += ugp[densityIdx(nmat, k)];
1257 : :
1258 : 2274000 : auto u = ugp[ncomp+velocityIdx(nmat,0)];
1259 : 2274000 : auto v = ugp[ncomp+velocityIdx(nmat,1)];
1260 : 2274000 : auto w = ugp[ncomp+velocityIdx(nmat,2)];
1261 : :
1262 : 2274000 : std::vector< tk::real > apk( nmat, 0.0 );
1263 [ + + ]: 6822000 : for (std::size_t k=0; k<nmat; ++k)
1264 : : {
1265 : 4548000 : apk[k] = ugp[ncomp+pressureIdx(nmat,k)];
1266 : 4548000 : p += apk[k];
1267 : : }
1268 : :
1269 [ + - ][ - - ]: 2274000 : std::vector< std::array< tk::real, 3 > > fl( ncomp );
1270 : :
1271 : : // conservative part of momentum flux
1272 : 2274000 : fl[momentumIdx(nmat, 0)][0] = ugp[momentumIdx(nmat, 0)] * u + p;
1273 : 2274000 : fl[momentumIdx(nmat, 1)][0] = ugp[momentumIdx(nmat, 1)] * u;
1274 : 2274000 : fl[momentumIdx(nmat, 2)][0] = ugp[momentumIdx(nmat, 2)] * u;
1275 : :
1276 : 2274000 : fl[momentumIdx(nmat, 0)][1] = ugp[momentumIdx(nmat, 0)] * v;
1277 : 2274000 : fl[momentumIdx(nmat, 1)][1] = ugp[momentumIdx(nmat, 1)] * v + p;
1278 : 2274000 : fl[momentumIdx(nmat, 2)][1] = ugp[momentumIdx(nmat, 2)] * v;
1279 : :
1280 : 2274000 : fl[momentumIdx(nmat, 0)][2] = ugp[momentumIdx(nmat, 0)] * w;
1281 : 2274000 : fl[momentumIdx(nmat, 1)][2] = ugp[momentumIdx(nmat, 1)] * w;
1282 : 2274000 : fl[momentumIdx(nmat, 2)][2] = ugp[momentumIdx(nmat, 2)] * w + p;
1283 : :
1284 [ + + ]: 6822000 : for (std::size_t k=0; k<nmat; ++k)
1285 : : {
1286 : : // conservative part of volume-fraction flux
1287 : 4548000 : fl[volfracIdx(nmat, k)][0] = 0.0;
1288 : 4548000 : fl[volfracIdx(nmat, k)][1] = 0.0;
1289 : 4548000 : fl[volfracIdx(nmat, k)][2] = 0.0;
1290 : :
1291 : : // conservative part of material continuity flux
1292 : 4548000 : fl[densityIdx(nmat, k)][0] = u * ugp[densityIdx(nmat, k)];
1293 : 4548000 : fl[densityIdx(nmat, k)][1] = v * ugp[densityIdx(nmat, k)];
1294 : 4548000 : fl[densityIdx(nmat, k)][2] = w * ugp[densityIdx(nmat, k)];
1295 : :
1296 : : // conservative part of material total-energy flux
1297 : 4548000 : auto hmat = ugp[energyIdx(nmat, k)] + apk[k];
1298 : 4548000 : fl[energyIdx(nmat, k)][0] = u * hmat;
1299 : 4548000 : fl[energyIdx(nmat, k)][1] = v * hmat;
1300 : 4548000 : fl[energyIdx(nmat, k)][2] = w * hmat;
1301 : : }
1302 : :
1303 : : // NEED TO RETURN m_ncomp flux vectors in fl, not 5
1304 : :
1305 : 2274000 : return fl;
1306 : : }
1307 : :
1308 : : //! \brief Boundary state function providing the left and right state of a
1309 : : //! face at Dirichlet boundaries
1310 : : //! \param[in] system Equation system index
1311 : : //! \param[in] ncomp Number of scalar components in this PDE system
1312 : : //! \param[in] ul Left (domain-internal) state
1313 : : //! \param[in] x X-coordinate at which to compute the states
1314 : : //! \param[in] y Y-coordinate at which to compute the states
1315 : : //! \param[in] z Z-coordinate at which to compute the states
1316 : : //! \param[in] t Physical time
1317 : : //! \param[in] fn Unit face normal
1318 : : //! \return Left and right states for all scalar components in this PDE
1319 : : //! system
1320 : : //! \note The function signature must follow tk::StateFn. For multimat, the
1321 : : //! left or right state is the vector of conserved quantities, followed by
1322 : : //! the vector of primitive quantities appended to it.
1323 : : static tk::StateFn::result_type
1324 : 48000 : dirichlet( ncomp_t system, ncomp_t ncomp, const std::vector< tk::real >& ul,
1325 : : tk::real x, tk::real y, tk::real z, tk::real t,
1326 : : const std::array< tk::real, 3 >& fn )
1327 : : {
1328 : 48000 : const auto nmat =
1329 : : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[system];
1330 : :
1331 : 48000 : auto ur = Problem::initialize( system, ncomp, x, y, z, t );
1332 : : Assert( ur.size() == ncomp, "Incorrect size for boundary state vector" );
1333 : :
1334 [ + - ]: 48000 : ur.resize(ul.size());
1335 : :
1336 : : tk::real rho(0.0);
1337 [ + + ]: 192000 : for (std::size_t k=0; k<nmat; ++k)
1338 : 144000 : rho += ur[densityIdx(nmat, k)];
1339 : :
1340 : : // get primitives in boundary state
1341 : :
1342 : : // velocity
1343 : 48000 : ur[ncomp+velocityIdx(nmat, 0)] = ur[momentumIdx(nmat, 0)] / rho;
1344 : 48000 : ur[ncomp+velocityIdx(nmat, 1)] = ur[momentumIdx(nmat, 1)] / rho;
1345 : 48000 : ur[ncomp+velocityIdx(nmat, 2)] = ur[momentumIdx(nmat, 2)] / rho;
1346 : :
1347 : : // determine the speed of sound of the majority material
1348 : : auto almax(0.0);
1349 : : std::size_t kmax(0);
1350 [ + + ]: 192000 : for (std::size_t k=0; k<nmat; ++k)
1351 : : {
1352 [ + + ]: 144000 : if (ul[volfracIdx(nmat, k)] > almax)
1353 : : {
1354 : : almax = ul[volfracIdx(nmat, k)];
1355 : : kmax = k;
1356 : : }
1357 : : }
1358 : :
1359 [ + - ]: 48000 : auto vn = tk::dot({{ul[ncomp+velocityIdx(nmat, 0)],
1360 : : ul[ncomp+velocityIdx(nmat, 1)], ul[ncomp+velocityIdx(nmat, 2)]}}, fn);
1361 [ + - ]: 48000 : auto Ml = vn/eos_soundspeed< tag::multimat >(system,
1362 : : ul[densityIdx(nmat,kmax)], ul[ncomp+pressureIdx(nmat,kmax)],
1363 : : almax, kmax);
1364 : :
1365 : : // material pressures
1366 [ - + ]: 48000 : if (Ml > 1.0)
1367 : : {
1368 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
1369 : : {
1370 [ - - ]: 0 : tk::real arhomat = ur[densityIdx(nmat, k)];
1371 : 0 : tk::real arhoemat = ur[energyIdx(nmat, k)];
1372 : 0 : tk::real alphamat = ur[volfracIdx(nmat, k)];
1373 [ - - ]: 0 : ur[ncomp+pressureIdx(nmat, k)] = eos_pressure< tag::multimat >( system,
1374 : : arhomat, ur[ncomp+velocityIdx(nmat, 0)],
1375 : : ur[ncomp+velocityIdx(nmat, 1)], ur[ncomp+velocityIdx(nmat, 2)],
1376 : : arhoemat, alphamat, k );
1377 : : }
1378 : : }
1379 : : else
1380 : : {
1381 [ + + ]: 192000 : for (std::size_t k=0; k<nmat; ++k)
1382 : : {
1383 : 144000 : ur[ncomp+pressureIdx(nmat, k)] = ul[ncomp+pressureIdx(nmat, k)];
1384 : 144000 : ur[energyIdx(nmat, k)] = ur[volfracIdx(nmat, k)] *
1385 : 144000 : eos_totalenergy< tag::multimat >(system,
1386 : : ur[densityIdx(nmat, k)]/ur[volfracIdx(nmat, k)],
1387 : : ur[ncomp+velocityIdx(nmat, 0)],
1388 : : ur[ncomp+velocityIdx(nmat, 1)],
1389 : : ur[ncomp+velocityIdx(nmat, 2)],
1390 : : ul[ncomp+pressureIdx(nmat, k)]/ul[volfracIdx(nmat, k)], k);
1391 : : }
1392 : : }
1393 : :
1394 : : Assert( ur.size() == ncomp+nmat+3, "Incorrect size for appended "
1395 : : "boundary state vector" );
1396 : :
1397 [ + - ]: 48000 : return {{ std::move(ul), std::move(ur) }};
1398 : : }
1399 : :
1400 : : //! \brief Boundary state function providing the left and right state of a
1401 : : //! face at symmetry boundaries
1402 : : //! \param[in] system Equation system index
1403 : : //! \param[in] ncomp Number of scalar components in this PDE system
1404 : : //! \param[in] ul Left (domain-internal) state
1405 : : //! \param[in] fn Unit face normal
1406 : : //! \return Left and right states for all scalar components in this PDE
1407 : : //! system
1408 : : //! \note The function signature must follow tk::StateFn. For multimat, the
1409 : : //! left or right state is the vector of conserved quantities, followed by
1410 : : //! the vector of primitive quantities appended to it.
1411 : : static tk::StateFn::result_type
1412 : 1671000 : symmetry( ncomp_t system, ncomp_t ncomp, const std::vector< tk::real >& ul,
1413 : : tk::real, tk::real, tk::real, tk::real,
1414 : : const std::array< tk::real, 3 >& fn )
1415 : : {
1416 : 1671000 : const auto nmat =
1417 : : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[system];
1418 : :
1419 : : Assert( ul.size() == ncomp+nmat+3, "Incorrect size for appended internal "
1420 : : "state vector" );
1421 : :
1422 : : tk::real rho(0.0);
1423 [ + + ]: 5636400 : for (std::size_t k=0; k<nmat; ++k)
1424 : 3965400 : rho += ul[densityIdx(nmat, k)];
1425 : :
1426 : 1671000 : auto ur = ul;
1427 : :
1428 : : // Internal cell velocity components
1429 : 1671000 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
1430 : 1671000 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
1431 : 1671000 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
1432 : : // Normal component of velocity
1433 : 1671000 : auto vnl = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
1434 : : // Ghost state velocity components
1435 : 1671000 : auto v1r = v1l - 2.0*vnl*fn[0];
1436 : 1671000 : auto v2r = v2l - 2.0*vnl*fn[1];
1437 : 1671000 : auto v3r = v3l - 2.0*vnl*fn[2];
1438 : : // Boundary condition
1439 [ + + ]: 5636400 : for (std::size_t k=0; k<nmat; ++k)
1440 : : {
1441 : 3965400 : ur[volfracIdx(nmat, k)] = ul[volfracIdx(nmat, k)];
1442 : 3965400 : ur[densityIdx(nmat, k)] = ul[densityIdx(nmat, k)];
1443 : 3965400 : ur[energyIdx(nmat, k)] = ul[energyIdx(nmat, k)];
1444 : : }
1445 : 1671000 : ur[momentumIdx(nmat, 0)] = rho * v1r;
1446 : 1671000 : ur[momentumIdx(nmat, 1)] = rho * v2r;
1447 : 1671000 : ur[momentumIdx(nmat, 2)] = rho * v3r;
1448 : :
1449 : : // Internal cell primitive quantities using the separately reconstructed
1450 : : // primitive quantities. This is used to get ghost state for primitive
1451 : : // quantities
1452 : :
1453 : : // velocity
1454 : 1671000 : ur[ncomp+velocityIdx(nmat, 0)] = v1r;
1455 : 1671000 : ur[ncomp+velocityIdx(nmat, 1)] = v2r;
1456 : 1671000 : ur[ncomp+velocityIdx(nmat, 2)] = v3r;
1457 : : // material pressures
1458 [ + + ]: 5636400 : for (std::size_t k=0; k<nmat; ++k)
1459 : 3965400 : ur[ncomp+pressureIdx(nmat, k)] = ul[ncomp+pressureIdx(nmat, k)];
1460 : :
1461 : : Assert( ur.size() == ncomp+nmat+3, "Incorrect size for appended boundary "
1462 : : "state vector" );
1463 : :
1464 [ + - ]: 1671000 : return {{ std::move(ul), std::move(ur) }};
1465 : : }
1466 : :
1467 : : //! \brief Boundary state function providing the left and right state of a
1468 : : //! face at subsonic outlet boundaries
1469 : : //! \param[in] system Equation system index
1470 : : //! \param[in] ncomp Number of scalar components in this PDE system
1471 : : //! \param[in] ul Left (domain-internal) state
1472 : : //! \return Left and right states for all scalar components in this PDE
1473 : : //! system
1474 : : //! \details The subsonic outlet boudary calculation, implemented here, is
1475 : : //! based on the characteristic theory of hyperbolic systems. For subsonic
1476 : : //! outlet flow, there is 1 incoming characteristic per material.
1477 : : //! Therefore, we calculate the ghost cell state by taking material
1478 : : //! pressure from the outside and other quantities from the internal cell.
1479 : : //! \note The function signature must follow tk::StateFn
1480 : : static tk::StateFn::result_type
1481 : 0 : subsonicOutlet( ncomp_t system,
1482 : : ncomp_t ncomp,
1483 : : const std::vector< tk::real >& ul,
1484 : : tk::real, tk::real, tk::real, tk::real,
1485 : : const std::array< tk::real, 3 >& )
1486 : : {
1487 : 0 : const auto nmat =
1488 : : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[system];
1489 : :
1490 : 0 : auto fp =
1491 : : g_inputdeck.get< tag::param, eq, tag::farfield_pressure >()[ system ];
1492 : :
1493 : : Assert( ul.size() == ncomp+nmat+3, "Incorrect size for appended internal "
1494 : : "state vector" );
1495 : :
1496 : 0 : auto ur = ul;
1497 : :
1498 : : // Internal cell velocity components
1499 : 0 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
1500 : 0 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
1501 : 0 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
1502 : : // Boundary condition
1503 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
1504 : : {
1505 : 0 : ur[energyIdx(nmat, k)] = ul[volfracIdx(nmat, k)] * eos_totalenergy< eq >
1506 : 0 : (system, ur[densityIdx(nmat, k)]/ul[volfracIdx(nmat, k)], v1l, v2l,
1507 : : v3l, fp, k);
1508 : : }
1509 : :
1510 : : // Internal cell primitive quantities using the separately reconstructed
1511 : : // primitive quantities. This is used to get ghost state for primitive
1512 : : // quantities
1513 : :
1514 : : // velocity
1515 : 0 : ur[ncomp+velocityIdx(nmat, 0)] = v1l;
1516 : 0 : ur[ncomp+velocityIdx(nmat, 1)] = v2l;
1517 : 0 : ur[ncomp+velocityIdx(nmat, 2)] = v3l;
1518 : : // material pressures
1519 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
1520 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ul[volfracIdx(nmat, k)] * fp;
1521 : :
1522 : : Assert( ur.size() == ncomp+nmat+3, "Incorrect size for appended boundary "
1523 : : "state vector" );
1524 : :
1525 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
1526 : : }
1527 : :
1528 : : //! \brief Boundary state function providing the left and right state of a
1529 : : //! face at extrapolation boundaries
1530 : : //! \param[in] ul Left (domain-internal) state
1531 : : //! \return Left and right states for all scalar components in this PDE
1532 : : //! system
1533 : : //! \note The function signature must follow tk::StateFn. For multimat, the
1534 : : //! left or right state is the vector of conserved quantities, followed by
1535 : : //! the vector of primitive quantities appended to it.
1536 : : static tk::StateFn::result_type
1537 : 75000 : extrapolate( ncomp_t, ncomp_t, const std::vector< tk::real >& ul,
1538 : : tk::real, tk::real, tk::real, tk::real,
1539 : : const std::array< tk::real, 3 >& )
1540 : : {
1541 : 75000 : return {{ ul, ul }};
1542 : : }
1543 : : };
1544 : :
1545 : : } // dg::
1546 : :
1547 : : } // inciter::
1548 : :
1549 : : #endif // DGMultiMat_h
|
