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 (with velocity equilibrium) using
12 : : discontinuous Galerkin 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 "Integrate/SolidTerms.hpp"
40 : : #include "RiemannChoice.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 : : #include "MultiMat/BCFunctions.hpp"
48 : : #include "MultiMat/MiscMultiMatFns.hpp"
49 : :
50 : : namespace inciter {
51 : :
52 : : extern ctr::InputDeck g_inputdeck;
53 : :
54 : : namespace dg {
55 : :
56 : : //! \brief MultiMat used polymorphically with tk::DGPDE
57 : : //! \details The template arguments specify policies and are used to configure
58 : : //! the behavior of the class. The policies are:
59 : : //! - Physics - physics configuration, see PDE/MultiMat/Physics.h
60 : : //! - Problem - problem configuration, see PDE/MultiMat/Problem.h
61 : : //! \note The default physics is Euler, set in inciter::deck::check_multimat()
62 : : template< class Physics, class Problem >
63 : : class MultiMat {
64 : :
65 : : private:
66 : : using eq = tag::multimat;
67 : :
68 : : public:
69 : : //! Constructor
70 : 130 : explicit MultiMat() :
71 : 130 : m_ncomp( g_inputdeck.get< tag::ncomp >() ),
72 : : m_nprim(nprim()),
73 : : m_riemann( multimatRiemannSolver(
74 : 130 : g_inputdeck.get< tag::flux >() ) )
75 : : {
76 : : // associate boundary condition configurations with state functions
77 [ + - ][ + - ]: 910 : brigand::for_each< ctr::bclist::Keys >( ConfigBC( m_bc,
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ][ - - ]
78 : : { dirichlet
79 : : , symmetry
80 : : , invalidBC // Inlet BC not implemented
81 : : , invalidBC // Outlet BC not implemented
82 : : , farfieldOutlet
83 : : , extrapolate } ) );
84 : :
85 : : // EoS initialization
86 [ + - ]: 130 : initializeMaterialEoS( m_mat_blk );
87 : 130 : }
88 : :
89 : : //! Find the number of primitive quantities required for this PDE system
90 : : //! \return The number of primitive quantities required to be stored for
91 : : //! this PDE system
92 : 201 : std::size_t nprim() const
93 : : {
94 : 201 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
95 : : const auto& solidx = inciter::g_inputdeck.get<
96 : 201 : tag::matidxmap, tag::solidx >();
97 : :
98 : : // individual material pressures and three velocity components
99 : 201 : std::size_t np(nmat+3);
100 : :
101 [ + + ]: 687 : for (std::size_t k=0; k<nmat; ++k) {
102 [ + + ]: 486 : if (solidx[k] > 0) {
103 : : // individual material Cauchy stress tensor components
104 : 15 : np += 6;
105 : : }
106 : : }
107 : :
108 : 201 : 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
114 : : {
115 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
116 : 0 : return nmat;
117 : : }
118 : :
119 : : //! Assign number of DOFs per equation in the PDE system
120 : : //! \param[in,out] numEqDof Array storing number of Dofs for each PDE
121 : : //! equation
122 : 71 : void numEquationDofs(std::vector< std::size_t >& numEqDof) const
123 : : {
124 : : // all equation-dofs initialized to ndofs first
125 [ + + ]: 779 : for (std::size_t i=0; i<m_ncomp; ++i) {
126 : 708 : numEqDof.push_back(g_inputdeck.get< tag::ndof >());
127 : : }
128 : :
129 : : // volume fractions are P0Pm (ndof = 1) for multi-material simulations
130 : 71 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
131 [ + - ]: 71 : if(nmat > 1)
132 [ + + ]: 218 : for (std::size_t k=0; k<nmat; ++k)
133 : 147 : numEqDof[volfracIdx(nmat, k)] = 1;
134 : 71 : }
135 : :
136 : : //! Determine elements that lie inside the user-defined IC box
137 : : //! \param[in] geoElem Element geometry array
138 : : //! \param[in] nielem Number of internal elements
139 : : //! \param[in,out] inbox List of nodes at which box user ICs are set for
140 : : //! each IC box
141 : 71 : void IcBoxElems( const tk::Fields& geoElem,
142 : : std::size_t nielem,
143 : : std::vector< std::unordered_set< std::size_t > >& inbox ) const
144 : : {
145 : 71 : tk::BoxElems< eq >(geoElem, nielem, inbox);
146 : 71 : }
147 : :
148 : : //! Initalize the compressible flow equations, prepare for time integration
149 : : //! \param[in] L Block diagonal mass matrix
150 : : //! \param[in] inpoel Element-node connectivity
151 : : //! \param[in] coord Array of nodal coordinates
152 : : //! \param[in] inbox List of elements at which box user ICs are set for
153 : : //! each IC box
154 : : //! \param[in] elemblkid Element ids associated with mesh block ids where
155 : : //! user ICs are set
156 : : //! \param[in,out] unk Array of unknowns
157 : : //! \param[in] t Physical time
158 : : //! \param[in] nielem Number of internal elements
159 : 71 : void initialize( const tk::Fields& L,
160 : : const std::vector< std::size_t >& inpoel,
161 : : const tk::UnsMesh::Coords& coord,
162 : : const std::vector< std::unordered_set< std::size_t > >& inbox,
163 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
164 : : elemblkid,
165 : : tk::Fields& unk,
166 : : tk::real t,
167 : : const std::size_t nielem ) const
168 : : {
169 [ + - ][ + - ]: 71 : tk::initialize( m_ncomp, m_mat_blk, L, inpoel, coord,
170 : : Problem::initialize, unk, t, nielem );
171 : :
172 : 71 : const auto rdof = g_inputdeck.get< tag::rdof >();
173 : 71 : const auto& ic = g_inputdeck.get< tag::ic >();
174 : 71 : const auto& icbox = ic.get< tag::box >();
175 : 71 : const auto& icmbk = ic.get< tag::meshblock >();
176 : :
177 : 71 : const auto& bgpre = ic.get< tag::pressure >();
178 : 71 : const auto& bgtemp = ic.get< tag::temperature >();
179 : :
180 : : // Set initial conditions inside user-defined IC boxes and mesh blocks
181 [ + - ]: 142 : std::vector< tk::real > s(m_ncomp, 0.0);
182 [ + + ]: 26445 : for (std::size_t e=0; e<nielem; ++e) {
183 : : // inside user-defined box
184 [ + + ]: 26374 : if (!icbox.empty()) {
185 : 2412 : std::size_t bcnt = 0;
186 [ + + ]: 4824 : for (const auto& b : icbox) { // for all boxes
187 [ + - ][ + - ]: 2412 : if (inbox.size() > bcnt && inbox[bcnt].find(e) != inbox[bcnt].end())
[ + + ][ + + ]
188 : : {
189 [ + - ]: 2092 : std::vector< tk::real > box
190 : 1046 : { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
191 : 1046 : b.template get< tag::ymin >(), b.template get< tag::ymax >(),
192 : 1046 : b.template get< tag::zmin >(), b.template get< tag::zmax >() };
193 : 1046 : auto V_ex = (box[1]-box[0]) * (box[3]-box[2]) * (box[5]-box[4]);
194 [ + + ]: 23906 : for (std::size_t c=0; c<m_ncomp; ++c) {
195 : 22860 : auto mark = c*rdof;
196 [ + - ]: 22860 : s[c] = unk(e,mark);
197 : : // set high-order DOFs to zero
198 [ + + ]: 91440 : for (std::size_t i=1; i<rdof; ++i)
199 [ + - ]: 68580 : unk(e,mark+i) = 0.0;
200 : : }
201 [ + - ]: 1046 : initializeBox<ctr::boxList>( m_mat_blk, V_ex, t, b, bgpre,
202 : : bgtemp, s );
203 : : // store box-initialization in solution vector
204 [ + + ]: 23906 : for (std::size_t c=0; c<m_ncomp; ++c) {
205 : 22860 : auto mark = c*rdof;
206 [ + - ]: 22860 : unk(e,mark) = s[c];
207 : : }
208 : : }
209 : 2412 : ++bcnt;
210 : : }
211 : : }
212 : :
213 : : // inside user-specified mesh blocks
214 [ - + ]: 26374 : if (!icmbk.empty()) {
215 [ - - ]: 0 : for (const auto& b : icmbk) { // for all blocks
216 : 0 : auto blid = b.get< tag::blockid >();
217 : 0 : auto V_ex = b.get< tag::volume >();
218 [ - - ][ - - ]: 0 : if (elemblkid.find(blid) != elemblkid.end()) {
219 [ - - ]: 0 : const auto& elset = tk::cref_find(elemblkid, blid);
220 [ - - ][ - - ]: 0 : if (elset.find(e) != elset.end()) {
221 [ - - ]: 0 : initializeBox<ctr::meshblockList>( m_mat_blk, V_ex, t, b,
222 : : bgpre, bgtemp, s );
223 : : // store initialization in solution vector
224 [ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) {
225 : 0 : auto mark = c*rdof;
226 [ - - ]: 0 : unk(e,mark) = s[c];
227 : : }
228 : : }
229 : : }
230 : : }
231 : : }
232 : : }
233 : 71 : }
234 : :
235 : : //! Compute the left hand side block-diagonal mass matrix
236 : : //! \param[in] geoElem Element geometry array
237 : : //! \param[in,out] l Block diagonal mass matrix
238 : 71 : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const {
239 : 71 : const auto ndof = g_inputdeck.get< tag::ndof >();
240 : : // Unlike Compflow and Transport, there is a weak reconstruction about
241 : : // conservative variable after limiting function which will require the
242 : : // size of left hand side vector to be rdof
243 : 71 : tk::mass( m_ncomp, ndof, geoElem, l );
244 : 71 : }
245 : :
246 : : //! Update the interface cells to first order dofs
247 : : //! \param[in] unk Array of unknowns
248 : : //! \param[in] nielem Number of internal elements
249 : : // //! \param[in,out] ndofel Array of dofs
250 : : //! \details This function resets the high-order terms in interface cells.
251 : 5430 : void updateInterfaceCells( tk::Fields& unk,
252 : : std::size_t nielem,
253 : : std::vector< std::size_t >& /*ndofel*/ ) const
254 : : {
255 : 5430 : auto intsharp =
256 : 5430 : g_inputdeck.get< tag::multimat, tag::intsharp >();
257 : : // If this cell is not material interface, return this function
258 [ + + ]: 5430 : if(not intsharp) return;
259 : :
260 : 405 : auto rdof = g_inputdeck.get< tag::rdof >();
261 : 405 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
262 : : const auto& solidx = g_inputdeck.get<
263 : 405 : tag::matidxmap, tag::solidx >();
264 : :
265 [ + + ]: 254685 : for (std::size_t e=0; e<nielem; ++e) {
266 [ + - ]: 508560 : std::vector< std::size_t > matInt(nmat, 0);
267 [ + - ]: 508560 : std::vector< tk::real > alAvg(nmat, 0.0);
268 [ + + ]: 762840 : for (std::size_t k=0; k<nmat; ++k)
269 [ + - ]: 508560 : alAvg[k] = unk(e, volfracDofIdx(nmat,k,rdof,0));
270 [ + - ]: 254280 : auto intInd = interfaceIndicator(nmat, alAvg, matInt);
271 : :
272 : : // interface cells cannot be high-order
273 [ + + ]: 254280 : if (intInd) {
274 : : //ndofel[e] = 1;
275 [ + + ]: 38256 : for (std::size_t k=0; k<nmat; ++k) {
276 [ + - ]: 25504 : if (matInt[k]) {
277 [ + + ]: 102016 : for (std::size_t i=1; i<rdof; ++i) {
278 [ + - ]: 76512 : unk(e, densityDofIdx(nmat,k,rdof,i)) = 0.0;
279 [ + - ]: 76512 : unk(e, energyDofIdx(nmat,k,rdof,i)) = 0.0;
280 : : }
281 [ + + ]: 25504 : if (solidx[k] > 0) {
282 [ + + ]: 15392 : for (std::size_t i=0; i<3; ++i)
283 [ + + ]: 46176 : for (std::size_t j=0; j<3; ++j)
284 [ + + ]: 138528 : for (std::size_t idof=1; idof<rdof; ++idof) {
285 [ + - ]: 103896 : unk(e, deformDofIdx(nmat,solidx[k],i,j,rdof,idof)) = 0.0;
286 : : }
287 : : }
288 : : }
289 : : }
290 [ + + ]: 51008 : for (std::size_t idir=0; idir<3; ++idir) {
291 [ + + ]: 153024 : for (std::size_t i=1; i<rdof; ++i) {
292 [ + - ]: 114768 : unk(e, momentumDofIdx(nmat,idir,rdof,i)) = 0.0;
293 : : }
294 : : }
295 : : }
296 : : }
297 : : }
298 : :
299 : : //! Update the primitives for this PDE system
300 : : //! \param[in] unk Array of unknowns
301 : : //! \param[in] L The left hand side block-diagonal mass matrix
302 : : //! \param[in] geoElem Element geometry array
303 : : //! \param[in,out] prim Array of primitives
304 : : //! \param[in] nielem Number of internal elements
305 : : //! \details This function computes and stores the dofs for primitive
306 : : //! quantities, which are required for obtaining reconstructed states used
307 : : //! in the Riemann solver. See /PDE/Riemann/AUSM.hpp, where the
308 : : //! normal velocity for advection is calculated from independently
309 : : //! reconstructed velocities.
310 : 5501 : void updatePrimitives( const tk::Fields& unk,
311 : : const tk::Fields& L,
312 : : const tk::Fields& geoElem,
313 : : tk::Fields& prim,
314 : : std::size_t nielem ) const
315 : : {
316 : 5501 : const auto rdof = g_inputdeck.get< tag::rdof >();
317 : 5501 : const auto ndof = g_inputdeck.get< tag::ndof >();
318 : 5501 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
319 : 5501 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
320 : :
321 [ - + ][ - - ]: 5501 : Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
322 : : "vector and primitive vector at recent time step incorrect" );
323 [ - + ][ - - ]: 5501 : Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
324 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
325 [ - + ][ - - ]: 5501 : Assert( prim.nprop() == rdof*m_nprim, "Number of components in vector of "
[ - - ][ - - ]
[ - - ]
326 : : "primitive quantities must equal "+ std::to_string(rdof*m_nprim) );
327 : :
328 [ + + ]: 2379615 : for (std::size_t e=0; e<nielem; ++e)
329 : : {
330 [ + - ]: 4748228 : std::vector< tk::real > R(m_nprim*ndof, 0.0);
331 : :
332 [ + - ]: 2374114 : auto ng = tk::NGvol(ndof);
333 : :
334 : : // arrays for quadrature points
335 : 4748228 : std::array< std::vector< tk::real >, 3 > coordgp;
336 : 4748228 : std::vector< tk::real > wgp;
337 : :
338 [ + - ]: 2374114 : coordgp[0].resize( ng );
339 [ + - ]: 2374114 : coordgp[1].resize( ng );
340 [ + - ]: 2374114 : coordgp[2].resize( ng );
341 [ + - ]: 2374114 : wgp.resize( ng );
342 : :
343 [ + - ]: 2374114 : tk::GaussQuadratureTet( ng, coordgp, wgp );
344 : :
345 : : // Loop over quadrature points in element e
346 [ + + ]: 6779668 : for (std::size_t igp=0; igp<ng; ++igp)
347 : : {
348 : : // Compute the basis function
349 [ + - ]: 8811108 : auto B =
350 : 4405554 : tk::eval_basis( ndof, coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
351 : :
352 [ + - ]: 4405554 : auto w = wgp[igp] * geoElem(e, 0);
353 : :
354 [ + - ]: 8811108 : auto state = tk::eval_state( m_ncomp, rdof, ndof, e, unk, B );
355 : :
356 : : // bulk density at quadrature point
357 : 4405554 : tk::real rhob(0.0);
358 [ + + ]: 14316848 : for (std::size_t k=0; k<nmat; ++k)
359 : 9911294 : rhob += state[densityIdx(nmat, k)];
360 : :
361 : : // velocity vector at quadrature point
362 : : std::array< tk::real, 3 >
363 : 4405554 : vel{ state[momentumIdx(nmat, 0)]/rhob,
364 : 4405554 : state[momentumIdx(nmat, 1)]/rhob,
365 : 4405554 : state[momentumIdx(nmat, 2)]/rhob };
366 : :
367 [ + - ]: 8811108 : std::vector< tk::real > pri(m_nprim, 0.0);
368 : :
369 : : // Evaluate material pressure at quadrature point
370 [ + + ]: 14316848 : for(std::size_t imat = 0; imat < nmat; imat++)
371 : : {
372 : 9911294 : auto alphamat = state[volfracIdx(nmat, imat)];
373 : 9911294 : auto arhomat = state[densityIdx(nmat, imat)];
374 : 9911294 : auto arhoemat = state[energyIdx(nmat, imat)];
375 [ + - ]: 9911294 : auto gmat = getDeformGrad(nmat, imat, state);
376 : 9911294 : pri[pressureIdx(nmat,imat)] = m_mat_blk[imat].compute<
377 [ + - ]: 9911294 : EOS::pressure >( arhomat, vel[0], vel[1], vel[2], arhoemat,
378 : : alphamat, imat, gmat );
379 : :
380 [ + - ]: 9911294 : pri[pressureIdx(nmat,imat)] = constrain_pressure( m_mat_blk,
381 : : pri[pressureIdx(nmat,imat)], arhomat, alphamat, imat);
382 : :
383 [ + + ]: 9911294 : if (solidx[imat] > 0) {
384 [ + - ]: 102548 : auto asigmat = m_mat_blk[imat].computeTensor< EOS::CauchyStress >(
385 : : arhomat, vel[0], vel[1], vel[2], arhoemat,
386 : : alphamat, imat, gmat );
387 : :
388 : 102548 : pri[stressIdx(nmat,solidx[imat],0)] = asigmat[0][0];
389 : 102548 : pri[stressIdx(nmat,solidx[imat],1)] = asigmat[1][1];
390 : 102548 : pri[stressIdx(nmat,solidx[imat],2)] = asigmat[2][2];
391 : 102548 : pri[stressIdx(nmat,solidx[imat],3)] = asigmat[0][1];
392 : 102548 : pri[stressIdx(nmat,solidx[imat],4)] = asigmat[0][2];
393 : 102548 : pri[stressIdx(nmat,solidx[imat],5)] = asigmat[1][2];
394 : : }
395 : : }
396 : :
397 : : // Evaluate bulk velocity at quadrature point
398 [ + + ]: 17622216 : for (std::size_t idir=0; idir<3; ++idir) {
399 : 13216662 : pri[velocityIdx(nmat,idir)] = vel[idir];
400 : : }
401 : :
402 [ + + ]: 28148798 : for(std::size_t k = 0; k < m_nprim; k++)
403 : : {
404 : 23743244 : auto mark = k * ndof;
405 [ + + ]: 85575988 : for(std::size_t idof = 0; idof < ndof; idof++)
406 : 61832744 : R[mark+idof] += w * pri[k] * B[idof];
407 : : }
408 : : }
409 : :
410 : : // Update the DG solution of primitive variables
411 [ + + ]: 15960158 : for(std::size_t k = 0; k < m_nprim; k++)
412 : : {
413 : 13586044 : auto mark = k * ndof;
414 : 13586044 : auto rmark = k * rdof;
415 [ + + ]: 34789988 : for(std::size_t idof = 0; idof < ndof; idof++)
416 : : {
417 [ + - ][ + - ]: 21203944 : prim(e, rmark+idof) = R[mark+idof] / L(e, mark+idof);
418 [ + - ][ + + ]: 21203944 : if(fabs(prim(e, rmark+idof)) < 1e-16)
419 [ + - ]: 5745520 : prim(e, rmark+idof) = 0;
420 : : }
421 : : }
422 : : }
423 : 5501 : }
424 : :
425 : : //! Clean up the state of trace materials for this PDE system
426 : : //! \param[in] t Physical time
427 : : //! \param[in] geoElem Element geometry array
428 : : //! \param[in,out] unk Array of unknowns
429 : : //! \param[in,out] prim Array of primitives
430 : : //! \param[in] nielem Number of internal elements
431 : : //! \details This function cleans up the state of materials present in trace
432 : : //! quantities in each cell. Specifically, the state of materials with
433 : : //! very low volume-fractions in a cell is replaced by the state of the
434 : : //! material which is present in the largest quantity in that cell. This
435 : : //! becomes necessary when shocks pass through cells which contain a very
436 : : //! small amount of material. The state of that tiny material might
437 : : //! become unphysical and cause solution to diverge; thus requiring such
438 : : //! a "reset".
439 : 5430 : void cleanTraceMaterial( tk::real t,
440 : : const tk::Fields& geoElem,
441 : : tk::Fields& unk,
442 : : tk::Fields& prim,
443 : : std::size_t nielem ) const
444 : : {
445 : 5430 : [[maybe_unused]] const auto rdof = g_inputdeck.get< tag::rdof >();
446 : 5430 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
447 : :
448 [ - + ][ - - ]: 5430 : Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
449 : : "vector and primitive vector at recent time step incorrect" );
450 [ - + ][ - - ]: 5430 : Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
451 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
452 [ - + ][ - - ]: 5430 : Assert( prim.nprop() == rdof*m_nprim, "Number of components in vector of "
[ - - ][ - - ]
[ - - ]
453 : : "primitive quantities must equal "+ std::to_string(rdof*m_nprim) );
454 : :
455 : 5430 : auto neg_density = cleanTraceMultiMat(t, nielem, m_mat_blk, geoElem, nmat,
456 : : unk, prim);
457 : :
458 [ - + ][ - - ]: 5430 : if (neg_density) Throw("Negative partial density.");
[ - - ][ - - ]
459 : 5430 : }
460 : :
461 : : //! Reconstruct second-order solution from first-order
462 : : //! \param[in] geoElem Element geometry array
463 : : //! \param[in] fd Face connectivity and boundary conditions object
464 : : //! \param[in] esup Elements-surrounding-nodes connectivity
465 : : //! \param[in] inpoel Element-node connectivity
466 : : //! \param[in] coord Array of nodal coordinates
467 : : //! \param[in,out] U Solution vector at recent time step
468 : : //! \param[in,out] P Vector of primitives at recent time step
469 : 4305 : void reconstruct( tk::real,
470 : : const tk::Fields&,
471 : : const tk::Fields& geoElem,
472 : : const inciter::FaceData& fd,
473 : : const std::map< std::size_t, std::vector< std::size_t > >&
474 : : esup,
475 : : const std::vector< std::size_t >& inpoel,
476 : : const tk::UnsMesh::Coords& coord,
477 : : tk::Fields& U,
478 : : tk::Fields& P ) const
479 : : {
480 : 4305 : const auto rdof = g_inputdeck.get< tag::rdof >();
481 : 4305 : const auto ndof = g_inputdeck.get< tag::ndof >();
482 : :
483 : 4305 : bool is_p0p1(false);
484 [ + - ][ + + ]: 4305 : if (rdof == 4 && ndof == 1)
485 : 3525 : is_p0p1 = true;
486 : :
487 : 4305 : const auto nelem = fd.Esuel().size()/4;
488 : 4305 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
489 : :
490 [ - + ][ - - ]: 4305 : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
491 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
492 : :
493 : : //----- reconstruction of conserved quantities -----
494 : : //--------------------------------------------------
495 : : // specify how many variables need to be reconstructed
496 : 8610 : std::vector< std::size_t > vars;
497 [ + + ][ + - ]: 44670 : for (std::size_t c=0; c<m_ncomp; ++c) vars.push_back(c);
498 : : // If DG is applied, reconstruct only volume fractions
499 [ + + ][ + - ]: 4305 : if (!is_p0p1 && ndof > 1)
500 : : {
501 : 780 : vars.clear();
502 [ + + ][ + - ]: 2340 : for (std::size_t k=0; k<nmat; ++k) vars.push_back(volfracIdx(nmat, k));
503 : : }
504 : :
505 : : // 1. solve 3x3 least-squares system
506 [ + + ]: 1031745 : for (std::size_t e=0; e<nelem; ++e)
507 : : {
508 : : // Reconstruct second-order dofs of volume-fractions in Taylor space
509 : : // using nodal-stencils, for a good interface-normal estimate
510 [ + - ]: 1027440 : tk::recoLeastSqExtStencil( rdof, e, esup, inpoel, geoElem, U, vars );
511 : : }
512 : :
513 : : // 2. transform reconstructed derivatives to Dubiner dofs
514 [ + - ]: 4305 : tk::transform_P0P1(rdof, nelem, inpoel, coord, U, vars);
515 : :
516 : : //----- reconstruction of primitive quantities -----
517 : : //--------------------------------------------------
518 : : // For multimat, conserved and primitive quantities are reconstructed
519 : : // separately.
520 [ + + ]: 4305 : if (is_p0p1) {
521 : 3525 : vars.clear();
522 [ + + ][ + - ]: 22230 : for (std::size_t c=0; c<m_nprim; ++c) vars.push_back(c);
523 : :
524 : : // 1.
525 [ + + ]: 530685 : for (std::size_t e=0; e<nelem; ++e)
526 : : {
527 : : // Reconstruct second-order dofs of volume-fractions in Taylor space
528 : : // using nodal-stencils, for a good interface-normal estimate
529 [ + - ]: 527160 : tk::recoLeastSqExtStencil( rdof, e, esup, inpoel, geoElem, P, vars );
530 : : }
531 : :
532 : : // 2.
533 [ + - ]: 3525 : tk::transform_P0P1(rdof, nelem, inpoel, coord, P, vars);
534 : : }
535 : 4305 : }
536 : :
537 : : //! Limit second-order solution, and primitive quantities separately
538 : : //! \param[in] t Physical time
539 : : //! \param[in] geoFace Face geometry array
540 : : //! \param[in] geoElem Element geometry array
541 : : //! \param[in] fd Face connectivity and boundary conditions object
542 : : //! \param[in] esup Elements-surrounding-nodes connectivity
543 : : //! \param[in] inpoel Element-node connectivity
544 : : //! \param[in] coord Array of nodal coordinates
545 : : //! \param[in] ndofel Vector of local number of degrees of freedome
546 : : //! \param[in] gid Local->global node id map
547 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
548 : : //! global node ids (key)
549 : : //! \param[in] uNodalExtrm Chare-boundary nodal extrema for conservative
550 : : //! variables
551 : : //! \param[in] pNodalExtrm Chare-boundary nodal extrema for primitive
552 : : //! variables
553 : : //! \param[in] mtInv Inverse of Taylor mass matrix
554 : : //! \param[in,out] U Solution vector at recent time step
555 : : //! \param[in,out] P Vector of primitives at recent time step
556 : : //! \param[in,out] shockmarker Vector of shock-marker values
557 : 4305 : void limit( [[maybe_unused]] tk::real t,
558 : : const tk::Fields& geoFace,
559 : : const tk::Fields& geoElem,
560 : : const inciter::FaceData& fd,
561 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
562 : : const std::vector< std::size_t >& inpoel,
563 : : const tk::UnsMesh::Coords& coord,
564 : : const std::vector< std::size_t >& ndofel,
565 : : const std::vector< std::size_t >& gid,
566 : : const std::unordered_map< std::size_t, std::size_t >& bid,
567 : : const std::vector< std::vector<tk::real> >& uNodalExtrm,
568 : : const std::vector< std::vector<tk::real> >& pNodalExtrm,
569 : : const std::vector< std::vector<tk::real> >& mtInv,
570 : : tk::Fields& U,
571 : : tk::Fields& P,
572 : : std::vector< std::size_t >& shockmarker ) const
573 : : {
574 [ - + ][ - - ]: 4305 : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
575 : : "vector and primitive vector at recent time step incorrect" );
576 : :
577 : 4305 : const auto limiter = g_inputdeck.get< tag::limiter >();
578 : 4305 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
579 : 4305 : const auto rdof = g_inputdeck.get< tag::rdof >();
580 : : const auto& solidx = g_inputdeck.get<
581 : 4305 : tag::matidxmap, tag::solidx >();
582 : :
583 : : // limit vectors of conserved and primitive quantities
584 [ - + ]: 4305 : if (limiter == ctr::LimiterType::SUPERBEEP1)
585 : : {
586 : 0 : SuperbeeMultiMat_P1( fd.Esuel(), inpoel, ndofel,
587 : : coord, solidx, U, P, nmat );
588 : : }
589 [ + - ][ + - ]: 4305 : else if (limiter == ctr::LimiterType::VERTEXBASEDP1 && rdof == 4)
590 : : {
591 [ + - ]: 8610 : VertexBasedMultiMat_P1( esup, inpoel, ndofel, fd.Esuel().size()/4,
592 : 4305 : m_mat_blk, fd, geoFace, geoElem, coord, flux, solidx, U, P,
593 : : nmat, shockmarker );
594 : : }
595 [ - - ][ - - ]: 0 : else if (limiter == ctr::LimiterType::VERTEXBASEDP1 && rdof == 10)
596 : : {
597 [ - - ]: 0 : VertexBasedMultiMat_P2( esup, inpoel, ndofel, fd.Esuel().size()/4,
598 : 0 : m_mat_blk, fd, geoFace, geoElem, coord, gid, bid,
599 : : uNodalExtrm, pNodalExtrm, mtInv, flux, solidx, U, P, nmat,
600 : : shockmarker );
601 : : }
602 [ - - ]: 0 : else if (limiter != ctr::LimiterType::NOLIMITER)
603 : : {
604 [ - - ][ - - ]: 0 : Throw("Limiter type not configured for multimat.");
[ - - ]
605 : : }
606 : 4305 : }
607 : :
608 : : //! Apply CPL to the conservative variable solution for this PDE system
609 : : //! \param[in] prim Array of primitive variables
610 : : //! \param[in] geoElem Element geometry array
611 : : //! \param[in] inpoel Element-node connectivity
612 : : //! \param[in] coord Array of nodal coordinates
613 : : //! \param[in,out] unk Array of conservative variables
614 : : //! \param[in] nielem Number of internal elements
615 : : //! \details This function applies CPL to obtain consistent dofs for
616 : : //! conservative quantities based on the limited primitive quantities.
617 : : //! See Pandare et al. (2023). On the Design of Stable,
618 : : //! Consistent, and Conservative High-Order Methods for Multi-Material
619 : : //! Hydrodynamics. J Comp Phys, 112313.
620 : 180 : void CPL( const tk::Fields& prim,
621 : : const tk::Fields& geoElem,
622 : : const std::vector< std::size_t >& inpoel,
623 : : const tk::UnsMesh::Coords& coord,
624 : : tk::Fields& unk,
625 : : std::size_t nielem ) const
626 : : {
627 : 180 : [[maybe_unused]] const auto rdof = g_inputdeck.get< tag::rdof >();
628 : 180 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
629 : :
630 [ - + ][ - - ]: 180 : Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
631 : : "vector and primitive vector at recent time step incorrect" );
632 [ - + ][ - - ]: 180 : Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
633 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
634 [ - + ][ - - ]: 180 : Assert( prim.nprop() == rdof*m_nprim, "Number of components in vector of "
[ - - ][ - - ]
[ - - ]
635 : : "primitive quantities must equal "+ std::to_string(rdof*m_nprim) );
636 : :
637 : 180 : correctLimConservMultiMat(nielem, m_mat_blk, nmat, inpoel,
638 : : coord, geoElem, prim, unk);
639 : 180 : }
640 : :
641 : : //! Return cell-average deformation gradient tensor
642 : : //! \param[in] unk Solution vector at recent time step
643 : : //! \param[in] nielem Number of internal elements
644 : : //! \details This function returns the bulk cell-average inverse
645 : : //! deformation gradient tensor
646 : 171 : std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
647 : : const tk::Fields& unk,
648 : : std::size_t nielem ) const
649 : : {
650 : 171 : const auto rdof = g_inputdeck.get< tag::rdof >();
651 : 171 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
652 : 171 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
653 : :
654 : 171 : std::array< std::vector< tk::real >, 9 > gb;
655 [ + - ][ + + ]: 171 : if (inciter::haveSolid(nmat, solidx)) {
656 [ + + ]: 140 : for (auto& gij : gb)
657 [ + - ]: 126 : gij.resize(nielem, 0.0);
658 [ + + ]: 6354 : for (std::size_t e=0; e<nielem; ++e) {
659 [ + + ]: 19020 : for (std::size_t k=0; k<nmat; ++k) {
660 [ + + ]: 12680 : if (solidx[k] > 0) {
661 [ + + ]: 32528 : for (std::size_t i=0; i<3; ++i)
662 [ + + ]: 97584 : for (std::size_t j=0; j<3; ++j)
663 [ + - ]: 146376 : gb[3*i+j][e] += unk(e, volfracDofIdx(nmat,k,rdof,0)) *
664 [ + - ]: 73188 : unk(e,deformDofIdx(nmat,solidx[k],i,j,rdof,0));
665 : : }
666 : : }
667 : : }
668 : : }
669 : :
670 : 171 : return gb;
671 : : }
672 : :
673 : :
674 : : //! Compute right hand side
675 : : //! \param[in] t Physical time
676 : : //! \param[in] geoFace Face geometry array
677 : : //! \param[in] geoElem Element geometry array
678 : : //! \param[in] fd Face connectivity and boundary conditions object
679 : : //! \param[in] inpoel Element-node connectivity
680 : : //! \param[in] coord Array of nodal coordinates
681 : : //! \param[in] U Solution vector at recent time step
682 : : //! \param[in] P Primitive vector at recent time step
683 : : //! \param[in] ndofel Vector of local number of degrees of freedom
684 : : //! \param[in] dt Delta time
685 : : //! \param[in,out] R Right-hand side vector computed
686 : 5430 : void rhs( tk::real t,
687 : : const tk::Fields& geoFace,
688 : : const tk::Fields& geoElem,
689 : : const inciter::FaceData& fd,
690 : : const std::vector< std::size_t >& inpoel,
691 : : const std::vector< std::unordered_set< std::size_t > >&,
692 : : const tk::UnsMesh::Coords& coord,
693 : : const tk::Fields& U,
694 : : const tk::Fields& P,
695 : : const std::vector< std::size_t >& ndofel,
696 : : const tk::real dt,
697 : : tk::Fields& R ) const
698 : : {
699 : 5430 : const auto ndof = g_inputdeck.get< tag::ndof >();
700 : 5430 : const auto rdof = g_inputdeck.get< tag::rdof >();
701 : 5430 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
702 : 5430 : const auto intsharp =
703 : 5430 : g_inputdeck.get< tag::multimat, tag::intsharp >();
704 : : const auto& solidx = inciter::g_inputdeck.get<
705 : 5430 : tag::matidxmap, tag::solidx >();
706 [ + - ]: 5430 : auto nsld = numSolids(nmat, solidx);
707 : :
708 : 5430 : const auto nelem = fd.Esuel().size()/4;
709 : :
710 [ - + ][ - - ]: 5430 : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
711 : : "vector and primitive vector at recent time step incorrect" );
712 [ - + ][ - - ]: 5430 : Assert( U.nunk() == R.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
713 : : "vector and right-hand side at recent time step incorrect" );
714 [ - + ][ - - ]: 5430 : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
715 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
716 [ - + ][ - - ]: 5430 : Assert( P.nprop() == rdof*m_nprim, "Number of components in primitive "
[ - - ][ - - ]
[ - - ]
717 : : "vector must equal "+ std::to_string(rdof*m_nprim) );
718 [ - + ][ - - ]: 5430 : Assert( R.nprop() == ndof*m_ncomp, "Number of components in right-hand "
[ - - ][ - - ]
[ - - ]
719 : : "side vector must equal "+ std::to_string(ndof*m_ncomp) );
720 [ - + ][ - - ]: 5430 : Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
[ - - ][ - - ]
721 : : "Mismatch in inpofa size" );
722 : :
723 : : // set rhs to zero
724 [ + - ]: 5430 : R.fill(0.0);
725 : :
726 : : // Allocate space for Riemann derivatives used in non-conservative terms.
727 : : // The following Riemann derivatives are stored, in order:
728 : : // 1) 3*nmat terms: derivatives of partial pressure of each material,
729 : : // for the energy equations.
730 : : // 2) ndof terms: derivatives of Riemann velocity times the basis
731 : : // function, for the volume fraction equations.
732 : : // 3) nmat*3*3*9 terms: 3 derivatives of u_l*g_ij for each material, for
733 : : // the deformation gradient equations.
734 : : // 4) 3*nsld terms: 3 derivatives of \alpha \sigma_ij for each solid
735 : : // material, for the energy equations.
736 : : std::vector< std::vector< tk::real > >
737 [ + - ][ + - ]: 16290 : riemannDeriv(3*nmat+ndof+3*nsld, std::vector<tk::real>(U.nunk(),0.0));
738 : :
739 : : // vectors to store the data of riemann velocity used for reconstruction
740 : : // in volume fraction equation
741 [ + - ]: 10860 : std::vector< std::vector< tk::real > > vriem( U.nunk() );
742 [ + - ]: 10860 : std::vector< std::vector< tk::real > > riemannLoc( U.nunk() );
743 : :
744 : : // configure a no-op lambda for prescribed velocity
745 : 10345515 : auto velfn = []( ncomp_t, tk::real, tk::real, tk::real, tk::real ){
746 : 10345515 : return tk::VelFn::result_type(); };
747 : :
748 : : // compute internal surface flux integrals
749 [ + - ]: 10860 : tk::surfInt( nmat, m_mat_blk, t, ndof, rdof, inpoel, solidx,
750 [ + - ]: 5430 : coord, fd, geoFace, geoElem, m_riemann, velfn, U, P, ndofel,
751 : : dt, R, vriem, riemannLoc, riemannDeriv, intsharp );
752 : :
753 : : // compute optional source term
754 [ + - ][ + - ]: 5430 : tk::srcInt( m_mat_blk, t, ndof, fd.Esuel().size()/4, inpoel,
755 : : coord, geoElem, Problem::src, ndofel, R, nmat );
756 : :
757 [ + + ]: 5430 : if(ndof > 1)
758 : : // compute volume integrals
759 [ + - ][ + - ]: 780 : tk::volInt( nmat, t, m_mat_blk, ndof, rdof, nelem,
[ + - ]
760 : : inpoel, coord, geoElem, flux, velfn, U, P, ndofel, R,
761 : : intsharp );
762 : :
763 : : // compute boundary surface flux integrals
764 [ + + ]: 38010 : for (const auto& b : m_bc)
765 : 32580 : tk::bndSurfInt( nmat, m_mat_blk, ndof, rdof,
766 : 32580 : b.first, fd, geoFace, geoElem, inpoel, coord, t,
767 [ + - ][ + - ]: 32580 : m_riemann, velfn, b.second, U, P, ndofel, R, vriem,
768 : : riemannLoc, riemannDeriv, intsharp );
769 : :
770 [ - + ][ - - ]: 5430 : Assert( riemannDeriv.size() == 3*nmat+ndof+3*nsld, "Size of "
[ - - ][ - - ]
771 : : "Riemann derivative vector incorrect" );
772 : :
773 : : // get derivatives from riemannDeriv
774 [ + + ]: 48570 : for (std::size_t k=0; k<riemannDeriv.size(); ++k)
775 : : {
776 [ - + ][ - - ]: 43140 : Assert( riemannDeriv[k].size() == U.nunk(), "Riemann derivative vector "
[ - - ][ - - ]
777 : : "for non-conservative terms has incorrect size" );
778 [ + + ]: 30041580 : for (std::size_t e=0; e<U.nunk(); ++e)
779 [ + - ]: 29998440 : riemannDeriv[k][e] /= geoElem(e, 0);
780 : : }
781 : :
782 : : // compute volume integrals of non-conservative terms
783 [ + - ]: 5430 : tk::nonConservativeInt( nmat, m_mat_blk, ndof, rdof, nelem,
784 : : inpoel, coord, geoElem, U, P, riemannDeriv,
785 : : ndofel, R, intsharp );
786 : :
787 : : // Code below commented until details about the form of these terms in the
788 : : // \alpha_k g_k equations are sorted out.
789 : : // // Compute integrals for inverse deformation in solid materials
790 : : // if (inciter::haveSolid(nmat, solidx))
791 : : // tk::solidTermsVolInt( nmat, m_mat_blk, ndof, rdof, nelem,
792 : : // inpoel, coord, geoElem, U, P, ndofel, dt, R);
793 : :
794 : : // compute finite pressure relaxation terms
795 [ + + ]: 5430 : if (g_inputdeck.get< tag::multimat, tag::prelax >())
796 : : {
797 : 450 : const auto ct = g_inputdeck.get< tag::multimat,
798 : 450 : tag::prelax_timescale >();
799 [ + - ]: 450 : tk::pressureRelaxationInt( nmat, m_mat_blk, ndof,
800 : : rdof, nelem, inpoel, coord, geoElem, U, P,
801 : : ndofel, ct, R, intsharp );
802 : : }
803 : 5430 : }
804 : :
805 : : //! Evaluate the adaptive indicator and mark the ndof for each element
806 : : //! \param[in] nunk Number of unknowns
807 : : //! \param[in] coord Array of nodal coordinates
808 : : //! \param[in] inpoel Element-node connectivity
809 : : //! \param[in] fd Face connectivity and boundary conditions object
810 : : //! \param[in] unk Array of unknowns
811 : : //! \param[in] prim Array of primitive quantities
812 : : //! \param[in] indicator p-refinement indicator type
813 : : //! \param[in] ndof Number of degrees of freedom in the solution
814 : : //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
815 : : //! \param[in] tolref Tolerance for p-refinement
816 : : //! \param[in,out] ndofel Vector of local number of degrees of freedome
817 : 0 : void eval_ndof( std::size_t nunk,
818 : : [[maybe_unused]] const tk::UnsMesh::Coords& coord,
819 : : [[maybe_unused]] const std::vector< std::size_t >& inpoel,
820 : : const inciter::FaceData& fd,
821 : : const tk::Fields& unk,
822 : : const tk::Fields& prim,
823 : : inciter::ctr::PrefIndicatorType indicator,
824 : : std::size_t ndof,
825 : : std::size_t ndofmax,
826 : : tk::real tolref,
827 : : std::vector< std::size_t >& ndofel ) const
828 : : {
829 : 0 : const auto& esuel = fd.Esuel();
830 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
831 : :
832 [ - - ]: 0 : if(indicator == inciter::ctr::PrefIndicatorType::SPECTRAL_DECAY)
833 : 0 : spectral_decay(nmat, nunk, esuel, unk, prim, ndof, ndofmax, tolref,
834 : : ndofel);
835 : : else
836 [ - - ][ - - ]: 0 : Throw( "No such adaptive indicator type" );
[ - - ]
837 : 0 : }
838 : :
839 : : //! Compute the minimum time step size
840 : : //! \param[in] fd Face connectivity and boundary conditions object
841 : : //! \param[in] geoFace Face geometry array
842 : : //! \param[in] geoElem Element geometry array
843 : : // //! \param[in] ndofel Vector of local number of degrees of freedom
844 : : //! \param[in] U Solution vector at recent time step
845 : : //! \param[in] P Vector of primitive quantities at recent time step
846 : : //! \param[in] nielem Number of internal elements
847 : : //! \return Minimum time step size
848 : : //! \details The allowable dt is calculated by looking at the maximum
849 : : //! wave-speed in elements surrounding each face, times the area of that
850 : : //! face. Once the maximum of this quantity over the mesh is determined,
851 : : //! the volume of each cell is divided by this quantity. A minimum of this
852 : : //! ratio is found over the entire mesh, which gives the allowable dt.
853 : 335 : tk::real dt( const std::array< std::vector< tk::real >, 3 >&,
854 : : const std::vector< std::size_t >&,
855 : : const inciter::FaceData& fd,
856 : : const tk::Fields& geoFace,
857 : : const tk::Fields& geoElem,
858 : : const std::vector< std::size_t >& /*ndofel*/,
859 : : const tk::Fields& U,
860 : : const tk::Fields& P,
861 : : const std::size_t nielem ) const
862 : : {
863 : 335 : const auto ndof = g_inputdeck.get< tag::ndof >();
864 : 335 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
865 : :
866 : 335 : auto mindt = timeStepSizeMultiMat( m_mat_blk, fd.Esuf(), geoFace, geoElem,
867 : : nielem, nmat, U, P);
868 : :
869 : 335 : tk::real dgp = 0.0;
870 [ + + ]: 335 : if (ndof == 4)
871 : : {
872 : 260 : dgp = 1.0;
873 : : }
874 [ - + ]: 75 : else if (ndof == 10)
875 : : {
876 : 0 : dgp = 2.0;
877 : : }
878 : :
879 : : // Scale smallest dt with CFL coefficient and the CFL is scaled by (2*p+1)
880 : : // where p is the order of the DG polynomial by linear stability theory.
881 : 335 : mindt /= (2.0*dgp + 1.0);
882 : 335 : return mindt;
883 : : }
884 : :
885 : : //! Extract the velocity field at cell nodes. Currently unused.
886 : : //! \param[in] U Solution vector at recent time step
887 : : //! \param[in] N Element node indices
888 : : //! \return Array of the four values of the velocity field
889 : : std::array< std::array< tk::real, 4 >, 3 >
890 : : velocity( const tk::Fields& U,
891 : : const std::array< std::vector< tk::real >, 3 >&,
892 : : const std::array< std::size_t, 4 >& N ) const
893 : : {
894 : : const auto rdof = g_inputdeck.get< tag::rdof >();
895 : : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
896 : :
897 : : std::array< std::array< tk::real, 4 >, 3 > v;
898 : : v[0] = U.extract( momentumDofIdx(nmat, 0, rdof, 0), N );
899 : : v[1] = U.extract( momentumDofIdx(nmat, 1, rdof, 0), N );
900 : : v[2] = U.extract( momentumDofIdx(nmat, 2, rdof, 0), N );
901 : :
902 : : std::vector< std::array< tk::real, 4 > > ar;
903 : : ar.resize(nmat);
904 : : for (std::size_t k=0; k<nmat; ++k)
905 : : ar[k] = U.extract( densityDofIdx(nmat, k, rdof, 0), N );
906 : :
907 : : std::array< tk::real, 4 > r{{ 0.0, 0.0, 0.0, 0.0 }};
908 : : for (std::size_t i=0; i<r.size(); ++i) {
909 : : for (std::size_t k=0; k<nmat; ++k)
910 : : r[i] += ar[k][i];
911 : : }
912 : :
913 : : std::transform( r.begin(), r.end(), v[0].begin(), v[0].begin(),
914 : : []( tk::real s, tk::real& d ){ return d /= s; } );
915 : : std::transform( r.begin(), r.end(), v[1].begin(), v[1].begin(),
916 : : []( tk::real s, tk::real& d ){ return d /= s; } );
917 : : std::transform( r.begin(), r.end(), v[2].begin(), v[2].begin(),
918 : : []( tk::real s, tk::real& d ){ return d /= s; } );
919 : : return v;
920 : : }
921 : :
922 : : //! Return a map that associates user-specified strings to functions
923 : : //! \return Map that associates user-specified strings to functions that
924 : : //! compute relevant quantities to be output to file
925 : 342 : std::map< std::string, tk::GetVarFn > OutVarFn() const
926 : 342 : { return MultiMatOutVarFn(); }
927 : :
928 : : //! Return analytic field names to be output to file
929 : : //! \return Vector of strings labelling analytic fields output in file
930 : 0 : std::vector< std::string > analyticFieldNames() const {
931 : 0 : auto nmat = g_inputdeck.get< eq, tag::nmat >();
932 : :
933 : 0 : return MultiMatFieldNames(nmat);
934 : : }
935 : :
936 : : //! Return time history field names to be output to file
937 : : //! \return Vector of strings labelling time history fields output in file
938 : 0 : std::vector< std::string > histNames() const {
939 : 0 : return MultiMatHistNames();
940 : : }
941 : :
942 : : //! Return surface field output going to file
943 : : std::vector< std::vector< tk::real > >
944 : 0 : surfOutput( const std::map< int, std::vector< std::size_t > >&,
945 : : tk::Fields& ) const
946 : : {
947 : 0 : std::vector< std::vector< tk::real > > s; // punt for now
948 : 0 : return s;
949 : : }
950 : :
951 : : //! Return time history field output evaluated at time history points
952 : : //! \param[in] h History point data
953 : : //! \param[in] inpoel Element-node connectivity
954 : : //! \param[in] coord Array of nodal coordinates
955 : : //! \param[in] U Array of unknowns
956 : : //! \param[in] P Array of primitive quantities
957 : : //! \return Vector of time history output of bulk flow quantities (density,
958 : : //! velocity, total energy, and pressure) evaluated at time history points
959 : : std::vector< std::vector< tk::real > >
960 : 0 : histOutput( const std::vector< HistData >& h,
961 : : const std::vector< std::size_t >& inpoel,
962 : : const tk::UnsMesh::Coords& coord,
963 : : const tk::Fields& U,
964 : : const tk::Fields& P ) const
965 : : {
966 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
967 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
968 : :
969 : 0 : const auto& x = coord[0];
970 : 0 : const auto& y = coord[1];
971 : 0 : const auto& z = coord[2];
972 : :
973 [ - - ]: 0 : std::vector< std::vector< tk::real > > Up(h.size());
974 : :
975 : 0 : std::size_t j = 0;
976 [ - - ]: 0 : for (const auto& p : h) {
977 : 0 : auto e = p.get< tag::elem >();
978 : 0 : auto chp = p.get< tag::coord >();
979 : :
980 : : // Evaluate inverse Jacobian
981 : 0 : std::array< std::array< tk::real, 3>, 4 > cp{{
982 : : {{ x[inpoel[4*e ]], y[inpoel[4*e ]], z[inpoel[4*e ]] }},
983 : 0 : {{ x[inpoel[4*e+1]], y[inpoel[4*e+1]], z[inpoel[4*e+1]] }},
984 : 0 : {{ x[inpoel[4*e+2]], y[inpoel[4*e+2]], z[inpoel[4*e+2]] }},
985 : 0 : {{ x[inpoel[4*e+3]], y[inpoel[4*e+3]], z[inpoel[4*e+3]] }} }};
986 : 0 : auto J = tk::inverseJacobian( cp[0], cp[1], cp[2], cp[3] );
987 : :
988 : : // evaluate solution at history-point
989 : 0 : std::array< tk::real, 3 > dc{{chp[0]-cp[0][0], chp[1]-cp[0][1],
990 : 0 : chp[2]-cp[0][2]}};
991 [ - - ]: 0 : auto B = tk::eval_basis(rdof, tk::dot(J[0],dc), tk::dot(J[1],dc),
992 : 0 : tk::dot(J[2],dc));
993 [ - - ]: 0 : auto uhp = eval_state(m_ncomp, rdof, rdof, e, U, B);
994 [ - - ]: 0 : auto php = eval_state(m_nprim, rdof, rdof, e, P, B);
995 : :
996 : : // store solution in history output vector
997 [ - - ]: 0 : Up[j].resize(6, 0.0);
998 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
999 : 0 : Up[j][0] += uhp[densityIdx(nmat,k)];
1000 : 0 : Up[j][4] += uhp[energyIdx(nmat,k)];
1001 : 0 : Up[j][5] += php[pressureIdx(nmat,k)];
1002 : : }
1003 : 0 : Up[j][1] = php[velocityIdx(nmat,0)];
1004 : 0 : Up[j][2] = php[velocityIdx(nmat,1)];
1005 : 0 : Up[j][3] = php[velocityIdx(nmat,2)];
1006 : 0 : ++j;
1007 : : }
1008 : :
1009 : 0 : return Up;
1010 : : }
1011 : :
1012 : : //! Return names of integral variables to be output to diagnostics file
1013 : : //! \return Vector of strings labelling integral variables output
1014 : 15 : std::vector< std::string > names() const
1015 : : {
1016 : 15 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
1017 : 15 : return MultiMatDiagNames(nmat);
1018 : : }
1019 : :
1020 : : //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
1021 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
1022 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
1023 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
1024 : : //! \param[in] t Physical time at which to evaluate the analytic solution
1025 : : //! \return Vector of analytic solution at given location and time
1026 : : std::vector< tk::real >
1027 : 0 : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
1028 : 0 : { return Problem::analyticSolution( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
1029 : :
1030 : : //! Return analytic solution for conserved variables
1031 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
1032 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
1033 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
1034 : : //! \param[in] t Physical time at which to evaluate the analytic solution
1035 : : //! \return Vector of analytic solution at given location and time
1036 : : std::vector< tk::real >
1037 : 1343008 : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
1038 : 1343008 : { return Problem::initialize( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
1039 : :
1040 : : //! Return cell-averaged specific total energy for an element
1041 : : //! \param[in] e Element id for which total energy is required
1042 : : //! \param[in] unk Vector of conserved quantities
1043 : : //! \return Cell-averaged specific total energy for given element
1044 : 842728 : tk::real sp_totalenergy(std::size_t e, const tk::Fields& unk) const
1045 : : {
1046 : 842728 : const auto rdof = g_inputdeck.get< tag::rdof >();
1047 : 842728 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
1048 : :
1049 : 842728 : tk::real sp_te(0.0);
1050 : : // sum each material total energy
1051 [ + + ]: 2918764 : for (std::size_t k=0; k<nmat; ++k) {
1052 : 2076036 : sp_te += unk(e, energyDofIdx(nmat,k,rdof,0));
1053 : : }
1054 : 842728 : return sp_te;
1055 : : }
1056 : :
1057 : : private:
1058 : : //! Number of components in this PDE system
1059 : : const ncomp_t m_ncomp;
1060 : : //! Number of primitive quantities stored in this PDE system
1061 : : const ncomp_t m_nprim;
1062 : : //! Riemann solver
1063 : : tk::RiemannFluxFn m_riemann;
1064 : : //! BC configuration
1065 : : BCStateFn m_bc;
1066 : : //! EOS material block
1067 : : std::vector< EOS > m_mat_blk;
1068 : :
1069 : : //! Evaluate conservative part of physical flux function for this PDE system
1070 : : //! \param[in] ncomp Number of scalar components in this PDE system
1071 : : //! \param[in] ugp Numerical solution at the Gauss point at which to
1072 : : //! evaluate the flux
1073 : : //! \return Flux vectors for all components in this PDE system
1074 : : //! \note The function signature must follow tk::FluxFn
1075 : : static tk::FluxFn::result_type
1076 : 7374000 : flux( ncomp_t ncomp,
1077 : : const std::vector< EOS >& mat_blk,
1078 : : const std::vector< tk::real >& ugp,
1079 : : const std::vector< std::array< tk::real, 3 > >& )
1080 : : {
1081 : 7374000 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
1082 : :
1083 : 7374000 : return tk::fluxTerms(ncomp, nmat, mat_blk, ugp);
1084 : : }
1085 : :
1086 : : //! \brief Boundary state function providing the left and right state of a
1087 : : //! face at Dirichlet boundaries
1088 : : //! \param[in] ncomp Number of scalar components in this PDE system
1089 : : //! \param[in] mat_blk EOS material block
1090 : : //! \param[in] ul Left (domain-internal) state
1091 : : //! \param[in] x X-coordinate at which to compute the states
1092 : : //! \param[in] y Y-coordinate at which to compute the states
1093 : : //! \param[in] z Z-coordinate at which to compute the states
1094 : : //! \param[in] t Physical time
1095 : : //! \return Left and right states for all scalar components in this PDE
1096 : : //! system
1097 : : //! \note The function signature must follow tk::StateFn. For multimat, the
1098 : : //! left or right state is the vector of conserved quantities, followed by
1099 : : //! the vector of primitive quantities appended to it.
1100 : : static tk::StateFn::result_type
1101 : 48000 : dirichlet( ncomp_t ncomp,
1102 : : const std::vector< EOS >& mat_blk,
1103 : : const std::vector< tk::real >& ul, tk::real x, tk::real y,
1104 : : tk::real z, tk::real t, const std::array< tk::real, 3 >& )
1105 : : {
1106 : 48000 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
1107 : : const auto& solidx = g_inputdeck.get<
1108 : 48000 : tag::matidxmap, tag::solidx >();
1109 : :
1110 [ + - ]: 48000 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
1111 : :
1112 [ + - ]: 96000 : auto ur = Problem::initialize( ncomp, mat_blk, x, y, z, t );
1113 [ - + ][ - - ]: 48000 : Assert( ur.size() == ncomp, "Incorrect size for boundary state vector" );
[ - - ][ - - ]
1114 : :
1115 [ + - ]: 48000 : ur.resize(ul.size());
1116 : :
1117 : 48000 : tk::real rho(0.0);
1118 [ + + ]: 192000 : for (std::size_t k=0; k<nmat; ++k)
1119 : 144000 : rho += ur[densityIdx(nmat, k)];
1120 : :
1121 : : // get primitives in boundary state
1122 : :
1123 : : // velocity
1124 : 48000 : ur[ncomp+velocityIdx(nmat, 0)] = ur[momentumIdx(nmat, 0)] / rho;
1125 : 48000 : ur[ncomp+velocityIdx(nmat, 1)] = ur[momentumIdx(nmat, 1)] / rho;
1126 : 48000 : ur[ncomp+velocityIdx(nmat, 2)] = ur[momentumIdx(nmat, 2)] / rho;
1127 : :
1128 : : // material pressures
1129 [ + + ]: 192000 : for (std::size_t k=0; k<nmat; ++k)
1130 : : {
1131 [ + - ]: 144000 : auto gk = getDeformGrad(nmat, k, ur);
1132 [ + - ]: 288000 : ur[ncomp+pressureIdx(nmat, k)] = mat_blk[k].compute< EOS::pressure >(
1133 : 144000 : ur[densityIdx(nmat, k)], ur[ncomp+velocityIdx(nmat, 0)],
1134 : 144000 : ur[ncomp+velocityIdx(nmat, 1)], ur[ncomp+velocityIdx(nmat, 2)],
1135 : 144000 : ur[energyIdx(nmat, k)], ur[volfracIdx(nmat, k)], k, gk );
1136 : : }
1137 : :
1138 [ - + ][ - - ]: 48000 : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
1139 : : "boundary state vector" );
1140 : :
1141 [ + - ]: 96000 : return {{ std::move(ul), std::move(ur) }};
1142 : : }
1143 : :
1144 : : // Other boundary condition types that do not depend on "Problem" should be
1145 : : // added in BCFunctions.hpp
1146 : : };
1147 : :
1148 : : } // dg::
1149 : :
1150 : : } // inciter::
1151 : :
1152 : : #endif // DGMultiMat_h
|
