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 : : #include "EoS/GetMatProp.hpp"
50 : :
51 : : namespace inciter {
52 : :
53 : : extern ctr::InputDeck g_inputdeck;
54 : :
55 : : namespace dg {
56 : :
57 : : //! \brief MultiMat used polymorphically with tk::DGPDE
58 : : //! \details The template arguments specify policies and are used to configure
59 : : //! the behavior of the class. The policies are:
60 : : //! - Physics - physics configuration, see PDE/MultiMat/Physics.h
61 : : //! - Problem - problem configuration, see PDE/MultiMat/Problem.h
62 : : //! \note The default physics is Euler, set in inciter::deck::check_multimat()
63 : : template< class Physics, class Problem >
64 : : class MultiMat {
65 : :
66 : : private:
67 : : using eq = tag::multimat;
68 : :
69 : : public:
70 : : //! Constructor
71 : 118 : explicit MultiMat() :
72 : 118 : m_ncomp( g_inputdeck.get< tag::ncomp >() ),
73 : : m_nprim(nprim()),
74 : : m_riemann( multimatRiemannSolver(
75 : 118 : g_inputdeck.get< tag::flux >() ) )
76 : : {
77 : : // associate boundary condition configurations with state functions
78 [ + - ][ + - ]: 1534 : brigand::for_each< ctr::bclist::Keys >( ConfigBC( m_bc,
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + + ]
[ + + ][ - - ]
[ - - ]
79 : : // BC State functions
80 : : { dirichlet
81 : : , symmetry
82 : : , invalidBC // Outlet BC not implemented
83 : : , farfield
84 : : , extrapolate
85 : : , noslipwall },
86 : : // BC Gradient functions
87 : : { noOpGrad
88 : : , symmetryGrad
89 : : , noOpGrad
90 : : , noOpGrad
91 : : , noOpGrad
92 : : , noOpGrad }
93 : : ) );
94 : :
95 : : // Inlet BC has a different structure than above BCs, so it must be
96 : : // handled differently than with ConfigBC
97 [ + - ][ + - ]: 118 : ConfigInletBC(m_bc, inlet, zeroGrad);
[ + - ]
98 : :
99 : : // EoS initialization
100 [ + - ]: 118 : initializeMaterialEoS( m_mat_blk );
101 : 118 : }
102 : :
103 : : //! Find the number of primitive quantities required for this PDE system
104 : : //! \return The number of primitive quantities required to be stored for
105 : : //! this PDE system
106 : 183 : std::size_t nprim() const
107 : : {
108 : 183 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
109 : : const auto& solidx = inciter::g_inputdeck.get<
110 : 183 : tag::matidxmap, tag::solidx >();
111 : :
112 : : // individual material pressures and three velocity components
113 : 183 : std::size_t np(nmat+3);
114 : :
115 [ + + ]: 633 : for (std::size_t k=0; k<nmat; ++k) {
116 [ - + ]: 450 : if (solidx[k] > 0) {
117 : : // individual material Cauchy stress tensor components
118 : 0 : np += 6;
119 : : }
120 : : }
121 : :
122 : 183 : return np;
123 : : }
124 : :
125 : : //! Find the number of materials set up for this PDE system
126 : : //! \return The number of materials set up for this PDE system
127 : 0 : std::size_t nmat() const
128 : : {
129 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
130 : 0 : return nmat;
131 : : }
132 : :
133 : : //! Assign number of DOFs per equation in the PDE system
134 : : //! \param[in,out] numEqDof Array storing number of Dofs for each PDE
135 : : //! equation
136 : 65 : void numEquationDofs(std::vector< std::size_t >& numEqDof) const
137 : : {
138 : : // all equation-dofs initialized to ndofs first
139 [ + + ]: 665 : for (std::size_t i=0; i<m_ncomp; ++i) {
140 : 600 : numEqDof.push_back(g_inputdeck.get< tag::ndof >());
141 : : }
142 : :
143 : : // volume fractions are P0Pm (ndof = 1) for multi-material simulations
144 : 65 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
145 [ + - ]: 65 : if(nmat > 1)
146 [ + + ]: 200 : for (std::size_t k=0; k<nmat; ++k)
147 : 135 : numEqDof[volfracIdx(nmat, k)] = 1;
148 : 65 : }
149 : :
150 : : //! Determine elements that lie inside the user-defined IC box
151 : : //! \param[in] geoElem Element geometry array
152 : : //! \param[in] nielem Number of internal elements
153 : : //! \param[in,out] inbox List of nodes at which box user ICs are set for
154 : : //! each IC box
155 : 65 : void IcBoxElems( const tk::Fields& geoElem,
156 : : std::size_t nielem,
157 : : std::vector< std::unordered_set< std::size_t > >& inbox ) const
158 : : {
159 : 65 : tk::BoxElems< eq >(geoElem, nielem, inbox);
160 : 65 : }
161 : :
162 : : //! Find how many 'stiff equations', which are the inverse
163 : : //! deformation equations because of plasticity
164 : : //! \return number of stiff equations
165 : 390 : std::size_t nstiffeq() const
166 : : {
167 : 390 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
168 : 390 : std::size_t nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
169 : 390 : return 9*numSolids(nmat, solidx);
170 : : }
171 : :
172 : : //! Find how many 'non-stiff equations', which are the inverse
173 : : //! deformation equations because of plasticity
174 : : //! \return number of stiff equations
175 : 130 : std::size_t nnonstiffeq() const
176 : : {
177 : 130 : return m_ncomp-nstiffeq();
178 : : }
179 : :
180 : : //! Locate the stiff equations.
181 : : //! \param[out] stiffEqIdx list with pointers to stiff equations
182 : 0 : void setStiffEqIdx( std::vector< std::size_t >& stiffEqIdx ) const
183 : : {
184 [ - - ][ - - ]: 0 : stiffEqIdx.resize(nstiffeq(), 0);
185 : 0 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
186 : 0 : std::size_t nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
187 : 0 : std::size_t icnt = 0;
188 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
189 [ - - ]: 0 : if (solidx[k] > 0)
190 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
191 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
192 : : {
193 : 0 : stiffEqIdx[icnt] =
194 : 0 : inciter::deformIdx(nmat, solidx[k], i, j);
195 : 0 : icnt++;
196 : : }
197 : 0 : }
198 : :
199 : : //! Locate the nonstiff equations.
200 : : //! \param[out] nonStiffEqIdx list with pointers to nonstiff equations
201 : 0 : void setNonStiffEqIdx( std::vector< std::size_t >& nonStiffEqIdx ) const
202 : : {
203 [ - - ][ - - ]: 0 : nonStiffEqIdx.resize(nnonstiffeq(), 0);
204 [ - - ]: 0 : for (std::size_t icomp=0; icomp<nnonstiffeq(); icomp++)
205 : 0 : nonStiffEqIdx[icomp] = icomp;
206 : 0 : }
207 : :
208 : : //! Initialize the compressible flow equations, prepare for time integration
209 : : //! \param[in] L Block diagonal mass matrix
210 : : //! \param[in] inpoel Element-node connectivity
211 : : //! \param[in] coord Array of nodal coordinates
212 : : //! \param[in] inbox List of elements at which box user ICs are set for
213 : : //! each IC box
214 : : //! \param[in] elemblkid Element ids associated with mesh block ids where
215 : : //! user ICs are set
216 : : //! \param[in,out] unk Array of unknowns
217 : : //! \param[in] t Physical time
218 : : //! \param[in] nielem Number of internal elements
219 : 65 : void initialize( const tk::Fields& L,
220 : : const std::vector< std::size_t >& inpoel,
221 : : const tk::UnsMesh::Coords& coord,
222 : : const std::vector< std::unordered_set< std::size_t > >& inbox,
223 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
224 : : elemblkid,
225 : : tk::Fields& unk,
226 : : tk::real t,
227 : : const std::size_t nielem ) const
228 : : {
229 [ + - ][ + - ]: 65 : tk::initialize( m_ncomp, m_mat_blk, L, inpoel, coord,
230 : : Problem::initialize, unk, t, nielem );
231 : :
232 : 65 : const auto rdof = g_inputdeck.get< tag::rdof >();
233 : 65 : const auto& ic = g_inputdeck.get< tag::ic >();
234 : 65 : const auto& icbox = ic.get< tag::box >();
235 : 65 : const auto& icmbk = ic.get< tag::meshblock >();
236 : :
237 : 65 : const auto& bgpre = ic.get< tag::pressure >();
238 : 65 : const auto& bgtemp = ic.get< tag::temperature >();
239 : :
240 : : // Set initial conditions inside user-defined IC boxes and mesh blocks
241 [ + - ]: 130 : std::vector< tk::real > s(m_ncomp, 0.0);
242 [ + + ]: 22511 : for (std::size_t e=0; e<nielem; ++e) {
243 : : // inside user-defined box
244 [ - + ]: 22446 : if (!icbox.empty()) {
245 : 0 : std::size_t bcnt = 0;
246 [ - - ]: 0 : for (const auto& b : icbox) { // for all boxes
247 [ - - ][ - - ]: 0 : if (inbox.size() > bcnt && inbox[bcnt].find(e) != inbox[bcnt].end())
[ - - ][ - - ]
248 : : {
249 [ - - ]: 0 : std::vector< tk::real > box
250 : 0 : { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
251 : 0 : b.template get< tag::ymin >(), b.template get< tag::ymax >(),
252 : 0 : b.template get< tag::zmin >(), b.template get< tag::zmax >() };
253 : 0 : auto V_ex = (box[1]-box[0]) * (box[3]-box[2]) * (box[5]-box[4]);
254 [ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) {
255 : 0 : auto mark = c*rdof;
256 [ - - ]: 0 : s[c] = unk(e,mark);
257 : : // set high-order DOFs to zero
258 [ - - ]: 0 : for (std::size_t i=1; i<rdof; ++i)
259 [ - - ]: 0 : unk(e,mark+i) = 0.0;
260 : : }
261 [ - - ]: 0 : initializeBox<ctr::boxList>( m_mat_blk, V_ex, t, b, bgpre,
262 : : bgtemp, s );
263 : : // store box-initialization in solution vector
264 [ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) {
265 : 0 : auto mark = c*rdof;
266 [ - - ]: 0 : unk(e,mark) = s[c];
267 : : }
268 : : }
269 : 0 : ++bcnt;
270 : : }
271 : : }
272 : :
273 : : // inside user-specified mesh blocks
274 [ - + ]: 22446 : if (!icmbk.empty()) {
275 [ - - ]: 0 : for (const auto& b : icmbk) { // for all blocks
276 : 0 : auto blid = b.get< tag::blockid >();
277 : 0 : auto V_ex = b.get< tag::volume >();
278 [ - - ][ - - ]: 0 : if (elemblkid.find(blid) != elemblkid.end()) {
279 [ - - ]: 0 : const auto& elset = tk::cref_find(elemblkid, blid);
280 [ - - ][ - - ]: 0 : if (elset.find(e) != elset.end()) {
281 [ - - ]: 0 : initializeBox<ctr::meshblockList>( m_mat_blk, V_ex, t, b,
282 : : bgpre, bgtemp, s );
283 : : // store initialization in solution vector
284 [ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) {
285 : 0 : auto mark = c*rdof;
286 [ - - ]: 0 : unk(e,mark) = s[c];
287 : : }
288 : : }
289 : : }
290 : : }
291 : : }
292 : : }
293 : 65 : }
294 : :
295 : : //! Compute density constraint for a given material
296 : : //! \param[in] nelem Number of elements
297 : : //! \param[in] unk Array of unknowns
298 : : //! \param[out] densityConstr Density Constraint: rho/(rho0*det(g))
299 : 155 : void computeDensityConstr( std::size_t nelem,
300 : : tk::Fields& unk,
301 : : std::vector< tk::real >& densityConstr) const
302 : : {
303 : 155 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
304 : 155 : std::size_t rdof = g_inputdeck.get< tag::rdof >();
305 : 155 : std::size_t nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
306 [ + + ]: 64461 : for (std::size_t e=0; e<nelem; ++e)
307 : 64306 : densityConstr[e] = 0.0;
308 [ + + ]: 480 : for (std::size_t imat=0; imat<nmat; ++imat)
309 [ - + ]: 325 : if (solidx[imat] > 0)
310 : : {
311 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
312 : : {
313 : : // Retrieve unknowns
314 [ - - ]: 0 : tk::real arho = unk(e, densityDofIdx(nmat, imat, rdof, 0));
315 : : std::array< std::array< tk::real, 3 >, 3 > g;
316 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
317 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
318 [ - - ]: 0 : g[i][j] = unk(e, deformDofIdx(nmat, solidx[imat], i, j, rdof, 0));
319 : : // Compute determinant of g
320 : 0 : tk::real detg = tk::determinant(g);
321 : : // Compute constraint measure
322 [ - - ]: 0 : densityConstr[e] += arho/(m_mat_blk[imat].compute< EOS::rho0 >()*detg);
323 : : }
324 : : }
325 : : else
326 : : {
327 [ + + ]: 150795 : for (std::size_t e=0; e<nelem; ++e)
328 : : {
329 : : // Retrieve alpha and add it to the constraint measure
330 : 150470 : tk::real alpha = unk(e, volfracDofIdx(nmat, imat, rdof, 0));
331 : 150470 : densityConstr[e] += alpha;
332 : : }
333 : : }
334 : 155 : }
335 : :
336 : : //! Compute the left hand side block-diagonal mass matrix
337 : : //! \param[in] geoElem Element geometry array
338 : : //! \param[in,out] l Block diagonal mass matrix
339 : 65 : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const {
340 : 65 : const auto ndof = g_inputdeck.get< tag::ndof >();
341 : : // Unlike Compflow and Transport, there is a weak reconstruction about
342 : : // conservative variable after limiting function which will require the
343 : : // size of left hand side vector to be rdof
344 : 65 : tk::mass( m_ncomp, ndof, geoElem, l );
345 : 65 : }
346 : :
347 : : //! Update the interface cells to first order dofs
348 : : //! \param[in] unk Array of unknowns
349 : : //! \param[in] nielem Number of internal elements
350 : : //! \param[in,out] ndofel Array of dofs
351 : : //! \param[in,out] interface Vector of interface marker
352 : : //! \details This function resets the high-order terms in interface cells.
353 : 5205 : void updateInterfaceCells( tk::Fields& unk,
354 : : std::size_t nielem,
355 : : std::vector< std::size_t >& ndofel,
356 : : std::vector< std::size_t >& interface ) const
357 : : {
358 : 5205 : auto intsharp =
359 : 5205 : g_inputdeck.get< tag::multimat, tag::intsharp >();
360 : : // If this cell is not material interface, return this function
361 [ + + ]: 5205 : if(not intsharp) return;
362 : :
363 : 375 : auto rdof = g_inputdeck.get< tag::rdof >();
364 : 375 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
365 : : const auto& solidx = g_inputdeck.get<
366 : 375 : tag::matidxmap, tag::solidx >();
367 : :
368 [ + + ]: 227775 : for (std::size_t e=0; e<nielem; ++e) {
369 [ + - ]: 454800 : std::vector< std::size_t > matInt(nmat, 0);
370 [ + - ]: 454800 : std::vector< tk::real > alAvg(nmat, 0.0);
371 [ + + ]: 682200 : for (std::size_t k=0; k<nmat; ++k)
372 [ + - ]: 454800 : alAvg[k] = unk(e, volfracDofIdx(nmat,k,rdof,0));
373 [ + - ]: 227400 : auto intInd = interfaceIndicator(nmat, alAvg, matInt);
374 : :
375 : : // interface cells cannot be high-order
376 [ + + ]: 227400 : if (intInd) {
377 : 10834 : interface[e] = 1;
378 [ + + ]: 32502 : for (std::size_t k=0; k<nmat; ++k) {
379 [ + - ]: 21668 : if (matInt[k]) {
380 [ + + ]: 86672 : for (std::size_t i=1; i<rdof; ++i) {
381 [ + - ]: 65004 : unk(e, densityDofIdx(nmat,k,rdof,i)) = 0.0;
382 [ + - ]: 65004 : unk(e, energyDofIdx(nmat,k,rdof,i)) = 0.0;
383 : : }
384 [ - + ]: 21668 : if (solidx[k] > 0) {
385 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
386 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
387 [ - - ]: 0 : for (std::size_t idof=1; idof<rdof; ++idof) {
388 [ - - ]: 0 : unk(e, deformDofIdx(nmat,solidx[k],i,j,rdof,idof)) = 0.0;
389 : : }
390 : : }
391 : : }
392 : : }
393 [ + + ]: 43336 : for (std::size_t idir=0; idir<3; ++idir) {
394 [ + + ]: 130008 : for (std::size_t i=1; i<rdof; ++i) {
395 [ + - ]: 97506 : unk(e, momentumDofIdx(nmat,idir,rdof,i)) = 0.0;
396 : : }
397 : : }
398 : : } else {
399 : : // If the cell is marked as interface cell in the previous timestep
400 : : // and does not marked as interface for the current timestep, DGP2
401 : : // will be applied for the current timestep in p-adaptive process
402 : : // Please note this block is added since the spectral decay indicator
403 : : // does not applied to P0 cells.
404 [ + + ]: 216566 : if (interface[e] == 1) {
405 [ + - ][ - + ]: 6 : if(ndofel[e] < 10 && rdof == 10) {
[ - + ]
406 : 0 : ndofel[e] = 10;
407 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
408 [ - - ]: 0 : for (std::size_t i=1; i<rdof; ++i) {
409 [ - - ]: 0 : unk(e, densityDofIdx(nmat,k,rdof,i)) = 0.0;
410 [ - - ]: 0 : unk(e, energyDofIdx(nmat,k,rdof,i)) = 0.0;
411 : : }
412 : : }
413 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir) {
414 [ - - ]: 0 : for (std::size_t i=1; i<rdof; ++i) {
415 [ - - ]: 0 : unk(e, momentumDofIdx(nmat,idir,rdof,i)) = 0.0;
416 : : }
417 : : }
418 : : }
419 : : }
420 : 216566 : interface[e] = 0;
421 : : }
422 : : }
423 : : }
424 : :
425 : : //! Update the primitives for this PDE system
426 : : //! \param[in] unk Array of unknowns
427 : : //! \param[in] L The left hand side block-diagonal mass matrix
428 : : //! \param[in] geoElem Element geometry array
429 : : //! \param[in,out] prim Array of primitives
430 : : //! \param[in] nielem Number of internal elements
431 : : //! \param[in] ndofel Array of dofs
432 : : //! \details This function computes and stores the dofs for primitive
433 : : //! quantities, which are required for obtaining reconstructed states used
434 : : //! in the Riemann solver. See /PDE/Riemann/AUSM.hpp, where the
435 : : //! normal velocity for advection is calculated from independently
436 : : //! reconstructed velocities.
437 : 5270 : void updatePrimitives( const tk::Fields& unk,
438 : : const tk::Fields& L,
439 : : const tk::Fields& geoElem,
440 : : tk::Fields& prim,
441 : : std::size_t nielem,
442 : : std::vector< std::size_t >& ndofel ) const
443 : : {
444 : 5270 : const auto rdof = g_inputdeck.get< tag::rdof >();
445 : 5270 : const auto ndof = g_inputdeck.get< tag::ndof >();
446 : 5270 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
447 : 5270 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
448 : :
449 [ - + ][ - - ]: 5270 : Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
450 : : "vector and primitive vector at recent time step incorrect" );
451 [ - + ][ - - ]: 5270 : Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
452 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
453 [ - + ][ - - ]: 5270 : Assert( prim.nprop() == rdof*m_nprim, "Number of components in vector of "
[ - - ][ - - ]
[ - - ]
454 : : "primitive quantities must equal "+ std::to_string(rdof*m_nprim) );
455 : :
456 [ + + ]: 2189396 : for (std::size_t e=0; e<nielem; ++e)
457 : : {
458 [ + - ]: 4368252 : std::vector< tk::real > R(m_nprim*ndof, 0.0);
459 : :
460 [ + - ]: 2184126 : auto ng = tk::NGvol(ndof);
461 : :
462 : : // arrays for quadrature points
463 : 4368252 : std::array< std::vector< tk::real >, 3 > coordgp;
464 : 4368252 : std::vector< tk::real > wgp;
465 : :
466 [ + - ]: 2184126 : coordgp[0].resize( ng );
467 [ + - ]: 2184126 : coordgp[1].resize( ng );
468 [ + - ]: 2184126 : coordgp[2].resize( ng );
469 [ + - ]: 2184126 : wgp.resize( ng );
470 : :
471 [ + - ]: 2184126 : tk::GaussQuadratureTet( ng, coordgp, wgp );
472 : :
473 : : // Local degree of freedom
474 : 2184126 : auto dof_el = ndofel[e];
475 : :
476 : : // Loop over quadrature points in element e
477 [ + + ]: 6399692 : for (std::size_t igp=0; igp<ng; ++igp)
478 : : {
479 : : // Compute the basis function
480 [ + - ]: 8431132 : auto B =
481 : 4215566 : tk::eval_basis( dof_el, coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
482 : :
483 [ + - ]: 4215566 : auto w = wgp[igp] * geoElem(e, 0);
484 : :
485 [ + - ]: 8431132 : auto state = tk::eval_state( m_ncomp, rdof, dof_el, e, unk, B );
486 : :
487 : : // bulk density at quadrature point
488 : 4215566 : tk::real rhob(0.0);
489 [ + + ]: 13746884 : for (std::size_t k=0; k<nmat; ++k)
490 : 9531318 : rhob += state[densityIdx(nmat, k)];
491 : :
492 : : // velocity vector at quadrature point
493 : : std::array< tk::real, 3 >
494 : 4215566 : vel{ state[momentumIdx(nmat, 0)]/rhob,
495 : 4215566 : state[momentumIdx(nmat, 1)]/rhob,
496 : 4215566 : state[momentumIdx(nmat, 2)]/rhob };
497 : :
498 [ + - ]: 8431132 : std::vector< tk::real > pri(m_nprim, 0.0);
499 : :
500 : : // Evaluate material pressure at quadrature point
501 [ + + ]: 13746884 : for(std::size_t imat = 0; imat < nmat; imat++)
502 : : {
503 : 9531318 : auto alphamat = state[volfracIdx(nmat, imat)];
504 : 9531318 : auto arhomat = state[densityIdx(nmat, imat)];
505 : 9531318 : auto arhoemat = state[energyIdx(nmat, imat)];
506 [ + - ]: 9531318 : auto gmat = getDeformGrad(nmat, imat, state);
507 : 9531318 : pri[pressureIdx(nmat,imat)] = m_mat_blk[imat].compute<
508 [ + - ]: 9531318 : EOS::pressure >( arhomat, vel[0], vel[1], vel[2], arhoemat,
509 : : alphamat, imat, gmat );
510 : :
511 [ + - ]: 9531318 : pri[pressureIdx(nmat,imat)] = constrain_pressure( m_mat_blk,
512 : : pri[pressureIdx(nmat,imat)], arhomat, alphamat, imat);
513 : :
514 [ - + ]: 9531318 : if (solidx[imat] > 0) {
515 [ - - ]: 0 : auto asigmat = m_mat_blk[imat].computeTensor< EOS::CauchyStress >(
516 : : arhomat, vel[0], vel[1], vel[2], arhoemat,
517 : : alphamat, imat, gmat );
518 : :
519 : 0 : pri[stressIdx(nmat,solidx[imat],0)] = asigmat[0][0];
520 : 0 : pri[stressIdx(nmat,solidx[imat],1)] = asigmat[1][1];
521 : 0 : pri[stressIdx(nmat,solidx[imat],2)] = asigmat[2][2];
522 : 0 : pri[stressIdx(nmat,solidx[imat],3)] = asigmat[0][1];
523 : 0 : pri[stressIdx(nmat,solidx[imat],4)] = asigmat[0][2];
524 : 0 : pri[stressIdx(nmat,solidx[imat],5)] = asigmat[1][2];
525 : : }
526 : : }
527 : :
528 : : // Evaluate bulk velocity at quadrature point
529 [ + + ]: 16862264 : for (std::size_t idir=0; idir<3; ++idir) {
530 : 12646698 : pri[velocityIdx(nmat,idir)] = vel[idir];
531 : : }
532 : :
533 [ + + ]: 26393582 : for(std::size_t k = 0; k < m_nprim; k++)
534 : : {
535 : 22178016 : auto mark = k * ndof;
536 [ + + ]: 82445532 : for(std::size_t idof = 0; idof < dof_el; idof++)
537 : 60267516 : R[mark+idof] += w * pri[k] * B[idof];
538 : : }
539 : : }
540 : :
541 : : // Update the DG solution of primitive variables
542 [ + + ]: 14204942 : for(std::size_t k = 0; k < m_nprim; k++)
543 : : {
544 : 12020816 : auto mark = k * ndof;
545 : 12020816 : auto rmark = k * rdof;
546 [ + + ]: 31659532 : for(std::size_t idof = 0; idof < dof_el; idof++)
547 : : {
548 [ + - ][ + - ]: 19638716 : prim(e, rmark+idof) = R[mark+idof] / L(e, mark+idof);
549 [ + - ][ + + ]: 19638716 : if(fabs(prim(e, rmark+idof)) < 1e-16)
550 [ + - ]: 5168733 : prim(e, rmark+idof) = 0;
551 : : }
552 : : }
553 : : }
554 : 5270 : }
555 : :
556 : : //! Clean up the state of trace materials for this PDE system
557 : : //! \param[in] t Physical time
558 : : //! \param[in] geoElem Element geometry array
559 : : //! \param[in,out] unk Array of unknowns
560 : : //! \param[in,out] prim Array of primitives
561 : : //! \param[in] nielem Number of internal elements
562 : : //! \details This function cleans up the state of materials present in trace
563 : : //! quantities in each cell. Specifically, the state of materials with
564 : : //! very low volume-fractions in a cell is replaced by the state of the
565 : : //! material which is present in the largest quantity in that cell. This
566 : : //! becomes necessary when shocks pass through cells which contain a very
567 : : //! small amount of material. The state of that tiny material might
568 : : //! become unphysical and cause solution to diverge; thus requiring such
569 : : //! a "reset".
570 : 5205 : void cleanTraceMaterial( tk::real t,
571 : : const tk::Fields& geoElem,
572 : : tk::Fields& unk,
573 : : tk::Fields& prim,
574 : : std::size_t nielem ) const
575 : : {
576 : 5205 : [[maybe_unused]] const auto rdof = g_inputdeck.get< tag::rdof >();
577 : 5205 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
578 : :
579 [ - + ][ - - ]: 5205 : Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
580 : : "vector and primitive vector at recent time step incorrect" );
581 [ - + ][ - - ]: 5205 : Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
582 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
583 [ - + ][ - - ]: 5205 : Assert( prim.nprop() == rdof*m_nprim, "Number of components in vector of "
[ - - ][ - - ]
[ - - ]
584 : : "primitive quantities must equal "+ std::to_string(rdof*m_nprim) );
585 : :
586 : 5205 : auto neg_density = cleanTraceMultiMat(t, nielem, m_mat_blk, geoElem, nmat,
587 : : unk, prim);
588 : :
589 [ - + ][ - - ]: 5205 : if (neg_density) Throw("Negative partial density.");
[ - - ][ - - ]
590 : 5205 : }
591 : :
592 : : //! Reconstruct second-order solution from first-order
593 : : //! \param[in] geoElem Element geometry array
594 : : //! \param[in] fd Face connectivity and boundary conditions object
595 : : //! \param[in] esup Elements-surrounding-nodes connectivity
596 : : //! \param[in] inpoel Element-node connectivity
597 : : //! \param[in] coord Array of nodal coordinates
598 : : //! \param[in,out] U Solution vector at recent time step
599 : : //! \param[in,out] P Vector of primitives at recent time step
600 : : //! \param[in] pref Indicator for p-adaptive algorithm
601 : : //! \param[in] ndofel Vector of local number of degrees of freedome
602 : 4080 : void reconstruct( tk::real,
603 : : const tk::Fields&,
604 : : const tk::Fields& geoElem,
605 : : const inciter::FaceData& fd,
606 : : const std::map< std::size_t, std::vector< std::size_t > >&
607 : : esup,
608 : : const std::vector< std::size_t >& inpoel,
609 : : const tk::UnsMesh::Coords& coord,
610 : : tk::Fields& U,
611 : : tk::Fields& P,
612 : : const bool pref,
613 : : const std::vector< std::size_t >& ndofel ) const
614 : : {
615 : 4080 : const auto rdof = g_inputdeck.get< tag::rdof >();
616 : 4080 : const auto ndof = g_inputdeck.get< tag::ndof >();
617 : :
618 : 4080 : bool is_p0p1(false);
619 [ + - ][ + + ]: 4080 : if (rdof == 4 && ndof == 1)
620 : 3300 : is_p0p1 = true;
621 : :
622 : 4080 : const auto nelem = fd.Esuel().size()/4;
623 : 4080 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
624 : :
625 [ - + ][ - - ]: 4080 : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
626 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
627 : :
628 : : //----- reconstruction of conserved quantities -----
629 : : //--------------------------------------------------
630 : :
631 [ + + ]: 845460 : for (std::size_t e=0; e<nelem; ++e)
632 : : {
633 : : // 1. specify how many variables need to be reconstructed
634 : 1682760 : std::vector< std::size_t > vars;
635 : : // for p-adaptive DG
636 [ - + ]: 841380 : if (pref) {
637 : : // If DG is applied, reconstruct only volume fractions
638 [ - - ]: 0 : if(ndofel[e] > 1) {
639 [ - - ][ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) vars.push_back(volfracIdx(nmat, k));
640 : : }
641 : : else // If P0P1 is applied for this element
642 [ - - ][ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) vars.push_back(c);
643 : : }
644 : : else {
645 : : // for P0P1, reconstruct all variables
646 [ + + ]: 841380 : if (is_p0p1)
647 [ + + ][ + - ]: 3411000 : for (std::size_t c=0; c<m_ncomp; ++c) vars.push_back(c);
648 : : // for high-order DG, reconstruct only volume fractions
649 [ + - ]: 500280 : else if (ndof > 1)
650 [ + + ][ + - ]: 1500840 : for (std::size_t k=0; k<nmat; ++k) vars.push_back(volfracIdx(nmat, k));
651 : : }
652 : :
653 : : // 2. solve 3x3 least-squares system
654 : : // Reconstruct second-order dofs of volume-fractions in Taylor space
655 : : // using nodal-stencils, for a good interface-normal estimate
656 [ + - ]: 841380 : tk::recoLeastSqExtStencil( rdof, e, esup, inpoel, geoElem, U, vars );
657 : :
658 : : // 3. transform reconstructed derivatives to Dubiner dofs
659 [ + - ]: 841380 : tk::transform_P0P1( rdof, e, inpoel, coord, U, vars );
660 : : }
661 : :
662 : : //----- reconstruction of primitive quantities -----
663 : : //--------------------------------------------------
664 : : // For multimat, conserved and primitive quantities are reconstructed
665 : : // separately.
666 : :
667 [ + + ]: 845460 : for (std::size_t e=0; e<nelem; ++e)
668 : : {
669 : : // There are two conditions that requires the reconstruction of the
670 : : // primitive variables:
671 : : // 1. p-adaptive is triggered and P0P1 scheme is applied to specific
672 : : // elements
673 : : // 2. p-adaptive is not triggered and P0P1 scheme is applied to the
674 : : // whole computation domain
675 [ - + ][ - - ]: 841380 : if ((pref && ndofel[e] == 1) || (!pref && is_p0p1)) {
[ + - ][ + + ]
[ + + ]
676 : 682200 : std::vector< std::size_t > vars;
677 [ + + ][ + - ]: 2046600 : for (std::size_t c=0; c<m_nprim; ++c) vars.push_back(c);
678 : :
679 : : // 1.
680 : : // Reconstruct second-order dofs of volume-fractions in Taylor space
681 : : // using nodal-stencils, for a good interface-normal estimate
682 [ + - ]: 341100 : tk::recoLeastSqExtStencil( rdof, e, esup, inpoel, geoElem, P, vars );
683 : :
684 : : // 2.
685 [ + - ]: 341100 : tk::transform_P0P1(rdof, e, inpoel, coord, P, vars);
686 : : }
687 : : }
688 : :
689 : 4080 : }
690 : :
691 : : //! Limit second-order solution, and primitive quantities separately
692 : : //! \param[in] t Physical time
693 : : //! \param[in] pref Indicator for p-adaptive algorithm
694 : : //! \param[in] geoFace Face geometry array
695 : : //! \param[in] geoElem Element geometry array
696 : : //! \param[in] fd Face connectivity and boundary conditions object
697 : : //! \param[in] esup Elements-surrounding-nodes connectivity
698 : : //! \param[in] inpoel Element-node connectivity
699 : : //! \param[in] coord Array of nodal coordinates
700 : : //! \param[in] ndofel Vector of local number of degrees of freedome
701 : : //! \param[in] gid Local->global node id map
702 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
703 : : //! global node ids (key)
704 : : //! \param[in] uNodalExtrm Chare-boundary nodal extrema for conservative
705 : : //! variables
706 : : //! \param[in] pNodalExtrm Chare-boundary nodal extrema for primitive
707 : : //! variables
708 : : //! \param[in] mtInv Inverse of Taylor mass matrix
709 : : //! \param[in,out] U Solution vector at recent time step
710 : : //! \param[in,out] P Vector of primitives at recent time step
711 : : //! \param[in,out] shockmarker Vector of shock-marker values
712 : 4080 : void limit( [[maybe_unused]] tk::real t,
713 : : const bool pref,
714 : : const tk::Fields& geoFace,
715 : : const tk::Fields& geoElem,
716 : : const inciter::FaceData& fd,
717 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
718 : : const std::vector< std::size_t >& inpoel,
719 : : const tk::UnsMesh::Coords& coord,
720 : : const std::vector< std::size_t >& ndofel,
721 : : const std::vector< std::size_t >& gid,
722 : : const std::unordered_map< std::size_t, std::size_t >& bid,
723 : : const std::vector< std::vector<tk::real> >& uNodalExtrm,
724 : : const std::vector< std::vector<tk::real> >& pNodalExtrm,
725 : : const std::vector< std::vector<tk::real> >& mtInv,
726 : : tk::Fields& U,
727 : : tk::Fields& P,
728 : : std::vector< std::size_t >& shockmarker ) const
729 : : {
730 [ - + ][ - - ]: 4080 : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
731 : : "vector and primitive vector at recent time step incorrect" );
732 : :
733 : 4080 : const auto limiter = g_inputdeck.get< tag::limiter >();
734 : 4080 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
735 : 4080 : const auto rdof = g_inputdeck.get< tag::rdof >();
736 : : const auto& solidx = g_inputdeck.get<
737 : 4080 : tag::matidxmap, tag::solidx >();
738 : :
739 : : // limit vectors of conserved and primitive quantities
740 [ - + ]: 4080 : if (limiter == ctr::LimiterType::SUPERBEEP1)
741 : : {
742 : 0 : SuperbeeMultiMat_P1( fd.Esuel(), inpoel, ndofel,
743 : : coord, solidx, U, P, nmat );
744 : : }
745 [ + - ][ + - ]: 4080 : else if (limiter == ctr::LimiterType::VERTEXBASEDP1 && rdof == 4)
746 : : {
747 [ + - ]: 8160 : VertexBasedMultiMat_P1( esup, inpoel, ndofel, fd.Esuel().size()/4,
748 : 4080 : m_mat_blk, fd, geoFace, geoElem, coord, flux, solidx, U, P,
749 : : nmat, shockmarker );
750 : : }
751 [ - - ][ - - ]: 0 : else if (limiter == ctr::LimiterType::VERTEXBASEDP1 && rdof == 10)
752 : : {
753 [ - - ]: 0 : VertexBasedMultiMat_P2( pref, esup, inpoel, ndofel, fd.Esuel().size()/4,
754 : 0 : m_mat_blk, fd, geoFace, geoElem, coord, gid, bid,
755 : : uNodalExtrm, pNodalExtrm, mtInv, flux, solidx, U, P, nmat,
756 : : shockmarker );
757 : : }
758 [ - - ]: 0 : else if (limiter != ctr::LimiterType::NOLIMITER)
759 : : {
760 [ - - ][ - - ]: 0 : Throw("Limiter type not configured for multimat.");
[ - - ]
761 : : }
762 : 4080 : }
763 : :
764 : : //! Apply CPL to the conservative variable solution for this PDE system
765 : : //! \param[in] prim Array of primitive variables
766 : : //! \param[in] geoElem Element geometry array
767 : : //! \param[in] inpoel Element-node connectivity
768 : : //! \param[in] coord Array of nodal coordinates
769 : : //! \param[in,out] unk Array of conservative variables
770 : : //! \param[in] nielem Number of internal elements
771 : : //! \details This function applies CPL to obtain consistent dofs for
772 : : //! conservative quantities based on the limited primitive quantities.
773 : : //! See Pandare et al. (2023). On the Design of Stable,
774 : : //! Consistent, and Conservative High-Order Methods for Multi-Material
775 : : //! Hydrodynamics. J Comp Phys, 112313.
776 : 30 : void CPL( const tk::Fields& prim,
777 : : const tk::Fields& geoElem,
778 : : const std::vector< std::size_t >& inpoel,
779 : : const tk::UnsMesh::Coords& coord,
780 : : tk::Fields& unk,
781 : : std::size_t nielem ) const
782 : : {
783 : 30 : [[maybe_unused]] const auto rdof = g_inputdeck.get< tag::rdof >();
784 : 30 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
785 : :
786 [ - + ][ - - ]: 30 : Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
787 : : "vector and primitive vector at recent time step incorrect" );
788 [ - + ][ - - ]: 30 : Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
789 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
790 [ - + ][ - - ]: 30 : Assert( prim.nprop() == rdof*m_nprim, "Number of components in vector of "
[ - - ][ - - ]
[ - - ]
791 : : "primitive quantities must equal "+ std::to_string(rdof*m_nprim) );
792 : :
793 : 30 : correctLimConservMultiMat(nielem, m_mat_blk, nmat, inpoel,
794 : : coord, geoElem, prim, unk);
795 : 30 : }
796 : :
797 : : //! Return cell-average deformation gradient tensor
798 : : //! \param[in] unk Solution vector at recent time step
799 : : //! \param[in] nielem Number of internal elements
800 : : //! \details This function returns the bulk cell-average inverse
801 : : //! deformation gradient tensor
802 : 0 : std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
803 : : const tk::Fields& unk,
804 : : std::size_t nielem ) const
805 : : {
806 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
807 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
808 : 0 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
809 : :
810 : 0 : std::array< std::vector< tk::real >, 9 > gb;
811 [ - - ][ - - ]: 0 : if (inciter::haveSolid(nmat, solidx)) {
812 [ - - ]: 0 : for (auto& gij : gb)
813 [ - - ]: 0 : gij.resize(nielem, 0.0);
814 [ - - ]: 0 : for (std::size_t e=0; e<nielem; ++e) {
815 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
816 [ - - ]: 0 : if (solidx[k] > 0) {
817 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
818 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
819 [ - - ]: 0 : gb[3*i+j][e] += unk(e, volfracDofIdx(nmat,k,rdof,0)) *
820 [ - - ]: 0 : unk(e,deformDofIdx(nmat,solidx[k],i,j,rdof,0));
821 : : }
822 : : }
823 : : }
824 : : }
825 : :
826 : 0 : return gb;
827 : : }
828 : :
829 : :
830 : : //! Reset the high order solution for p-adaptive scheme
831 : : //! \param[in] fd Face connectivity and boundary conditions object
832 : : //! \param[in,out] unk Solution vector at recent time step
833 : : //! \param[in,out] prim Primitive vector at recent time step
834 : : //! \param[in] ndofel Vector of local number of degrees of freedome
835 : : //! \details This function reset the high order coefficient for p-adaptive
836 : : //! solution polynomials. Unlike compflow class, the high order of fv
837 : : //! solution will not be reset since p0p1 is the base scheme for
838 : : //! multi-material p-adaptive DG method.
839 : 0 : void resetAdapSol( const inciter::FaceData& fd,
840 : : tk::Fields& unk,
841 : : tk::Fields& prim,
842 : : const std::vector< std::size_t >& ndofel ) const
843 : : {
844 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
845 : 0 : const auto ncomp = unk.nprop() / rdof;
846 : 0 : const auto nprim = prim.nprop() / rdof;
847 : :
848 [ - - ]: 0 : for(std::size_t e = 0; e < fd.Esuel().size()/4; e++)
849 : : {
850 [ - - ]: 0 : if(ndofel[e] < 10)
851 : : {
852 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c)
853 : : {
854 : 0 : auto mark = c*rdof;
855 : 0 : unk(e, mark+4) = 0.0;
856 : 0 : unk(e, mark+5) = 0.0;
857 : 0 : unk(e, mark+6) = 0.0;
858 : 0 : unk(e, mark+7) = 0.0;
859 : 0 : unk(e, mark+8) = 0.0;
860 : 0 : unk(e, mark+9) = 0.0;
861 : : }
862 [ - - ]: 0 : for (std::size_t c=0; c<nprim; ++c)
863 : : {
864 : 0 : auto mark = c*rdof;
865 : 0 : prim(e, mark+4) = 0.0;
866 : 0 : prim(e, mark+5) = 0.0;
867 : 0 : prim(e, mark+6) = 0.0;
868 : 0 : prim(e, mark+7) = 0.0;
869 : 0 : prim(e, mark+8) = 0.0;
870 : 0 : prim(e, mark+9) = 0.0;
871 : : }
872 : : }
873 : : }
874 : 0 : }
875 : :
876 : : //! Compute right hand side
877 : : //! \param[in] t Physical time
878 : : //! \param[in] pref Indicator for p-adaptive algorithm
879 : : //! \param[in] geoFace Face geometry array
880 : : //! \param[in] geoElem Element geometry array
881 : : //! \param[in] fd Face connectivity and boundary conditions object
882 : : //! \param[in] inpoel Element-node connectivity
883 : : //! \param[in] coord Array of nodal coordinates
884 : : //! \param[in] U Solution vector at recent time step
885 : : //! \param[in] P Primitive vector at recent time step
886 : : //! \param[in] ndofel Vector of local number of degrees of freedom
887 : : //! \param[in] dt Delta time
888 : : //! \param[in,out] R Right-hand side vector computed
889 : 5205 : void rhs( tk::real t,
890 : : const bool pref,
891 : : const tk::Fields& geoFace,
892 : : const tk::Fields& geoElem,
893 : : const inciter::FaceData& fd,
894 : : const std::vector< std::size_t >& inpoel,
895 : : const std::vector< std::unordered_set< std::size_t > >&,
896 : : const tk::UnsMesh::Coords& coord,
897 : : const tk::Fields& U,
898 : : const tk::Fields& P,
899 : : const std::vector< std::size_t >& ndofel,
900 : : const tk::real dt,
901 : : tk::Fields& R ) const
902 : : {
903 : 5205 : const auto ndof = g_inputdeck.get< tag::ndof >();
904 : 5205 : const auto rdof = g_inputdeck.get< tag::rdof >();
905 : 5205 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
906 : 5205 : const auto intsharp =
907 : 5205 : g_inputdeck.get< tag::multimat, tag::intsharp >();
908 : : const auto& solidx = inciter::g_inputdeck.get<
909 : 5205 : tag::matidxmap, tag::solidx >();
910 [ + - ]: 5205 : auto nsld = numSolids(nmat, solidx);
911 : :
912 : 5205 : const auto nelem = fd.Esuel().size()/4;
913 : :
914 [ - + ][ - - ]: 5205 : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
915 : : "vector and primitive vector at recent time step incorrect" );
916 [ - + ][ - - ]: 5205 : Assert( U.nunk() == R.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
917 : : "vector and right-hand side at recent time step incorrect" );
918 [ - + ][ - - ]: 5205 : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
919 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
920 [ - + ][ - - ]: 5205 : Assert( P.nprop() == rdof*m_nprim, "Number of components in primitive "
[ - - ][ - - ]
[ - - ]
921 : : "vector must equal "+ std::to_string(rdof*m_nprim) );
922 [ - + ][ - - ]: 5205 : Assert( R.nprop() == ndof*m_ncomp, "Number of components in right-hand "
[ - - ][ - - ]
[ - - ]
923 : : "side vector must equal "+ std::to_string(ndof*m_ncomp) );
924 [ - + ][ - - ]: 5205 : Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
[ - - ][ - - ]
925 : : "Mismatch in inpofa size" );
926 : :
927 : : // set rhs to zero
928 [ + - ]: 5205 : R.fill(0.0);
929 : :
930 : : // Allocate space for Riemann derivatives used in non-conservative terms.
931 : : // The following Riemann derivatives are stored, in order:
932 : : // 1) 3*nmat terms: derivatives of partial pressure of each material,
933 : : // for the energy equations.
934 : : // 2) ndof terms: derivatives of Riemann velocity times the basis
935 : : // function, for the volume fraction equations.
936 : : // 3) nmat*3*3*9 terms: 3 derivatives of u_l*g_ij for each material, for
937 : : // the deformation gradient equations.
938 : : // 4) 3*nsld terms: 3 derivatives of \alpha \sigma_ij for each solid
939 : : // material, for the energy equations.
940 : : std::vector< std::vector< tk::real > >
941 [ + - ][ + - ]: 15615 : riemannDeriv(3*nmat+ndof+3*nsld, std::vector<tk::real>(U.nunk(),0.0));
942 : :
943 : : // configure a no-op lambda for prescribed velocity
944 : 9919425 : auto velfn = []( ncomp_t, tk::real, tk::real, tk::real, tk::real ){
945 : 9919425 : return tk::VelFn::result_type(); };
946 : :
947 : : // compute internal surface flux integrals
948 [ + - ]: 10410 : tk::surfInt( pref, nmat, m_mat_blk, t, ndof, rdof, inpoel, solidx,
949 [ + - ]: 5205 : coord, fd, geoFace, geoElem, m_riemann, velfn, U, P, ndofel,
950 : : dt, R, riemannDeriv, intsharp );
951 : :
952 : : // compute optional source term
953 [ + - ][ + - ]: 5205 : tk::srcInt( m_mat_blk, t, ndof, fd.Esuel().size()/4, inpoel,
954 : : coord, geoElem, Problem::src, ndofel, R, nmat );
955 : :
956 [ + + ]: 5205 : if(ndof > 1)
957 : : // compute volume integrals
958 [ + - ][ + - ]: 780 : tk::volInt( nmat, t, m_mat_blk, ndof, rdof, nelem,
[ + - ]
959 : : inpoel, coord, geoElem, flux, velfn, U, P, ndofel, R,
960 : : intsharp );
961 : :
962 : : // compute boundary surface flux integrals
963 [ + + ]: 36435 : for (const auto& b : m_bc)
964 [ + - ][ + - ]: 62460 : tk::bndSurfInt( pref, nmat, m_mat_blk, ndof, rdof,
965 : 31230 : std::get<0>(b), fd, geoFace, geoElem, inpoel, coord, t,
966 : 62460 : m_riemann, velfn, std::get<1>(b), U, P, ndofel, R,
967 : : riemannDeriv, intsharp );
968 : :
969 [ - + ][ - - ]: 5205 : Assert( riemannDeriv.size() == 3*nmat+ndof+3*nsld, "Size of "
[ - - ][ - - ]
970 : : "Riemann derivative vector incorrect" );
971 : :
972 : : // get derivatives from riemannDeriv
973 [ + + ]: 46230 : for (std::size_t k=0; k<riemannDeriv.size(); ++k)
974 : : {
975 [ - + ][ - - ]: 41025 : Assert( riemannDeriv[k].size() == U.nunk(), "Riemann derivative vector "
[ - - ][ - - ]
976 : : "for non-conservative terms has incorrect size" );
977 [ + + ]: 28165725 : for (std::size_t e=0; e<U.nunk(); ++e)
978 [ + - ]: 28124700 : riemannDeriv[k][e] /= geoElem(e, 0);
979 : : }
980 : :
981 : : // compute volume integrals of non-conservative terms
982 [ + - ]: 5205 : tk::nonConservativeInt( pref, nmat, m_mat_blk, ndof, rdof, nelem,
983 : : inpoel, coord, geoElem, U, P, riemannDeriv,
984 : : ndofel, R, intsharp );
985 : :
986 : : // Compute integrals for inverse deformation correction in solid materials
987 [ + - ][ - + ]: 5205 : if (inciter::haveSolid(nmat, solidx) &&
[ - - ][ - + ]
988 : 0 : g_inputdeck.get< tag::multimat, tag::rho0constraint >())
989 [ - - ]: 0 : tk::solidTermsVolInt( nmat, m_mat_blk, ndof, rdof, nelem,
990 : : inpoel, coord, geoElem, U, P, ndofel,
991 : : dt, R);
992 : :
993 : : // compute finite pressure relaxation terms
994 [ + + ]: 5205 : if (g_inputdeck.get< tag::multimat, tag::prelax >())
995 : : {
996 : 375 : const auto ct = g_inputdeck.get< tag::multimat,
997 : 375 : tag::prelax_timescale >();
998 [ + - ]: 375 : tk::pressureRelaxationInt( pref, nmat, m_mat_blk, ndof,
999 : : rdof, nelem, inpoel, coord, geoElem, U, P,
1000 : : ndofel, ct, R, intsharp );
1001 : : }
1002 : 5205 : }
1003 : :
1004 : : //! Evaluate the adaptive indicator and mark the ndof for each element
1005 : : //! \param[in] nunk Number of unknowns
1006 : : //! \param[in] coord Array of nodal coordinates
1007 : : //! \param[in] inpoel Element-node connectivity
1008 : : //! \param[in] fd Face connectivity and boundary conditions object
1009 : : //! \param[in] unk Array of unknowns
1010 : : //! \param[in] prim Array of primitive quantities
1011 : : //! \param[in] indicator p-refinement indicator type
1012 : : //! \param[in] ndof Number of degrees of freedom in the solution
1013 : : //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
1014 : : //! \param[in] tolref Tolerance for p-refinement
1015 : : //! \param[in,out] ndofel Vector of local number of degrees of freedome
1016 : 0 : void eval_ndof( std::size_t nunk,
1017 : : [[maybe_unused]] const tk::UnsMesh::Coords& coord,
1018 : : [[maybe_unused]] const std::vector< std::size_t >& inpoel,
1019 : : const inciter::FaceData& fd,
1020 : : const tk::Fields& unk,
1021 : : [[maybe_unused]] const tk::Fields& prim,
1022 : : inciter::ctr::PrefIndicatorType indicator,
1023 : : std::size_t ndof,
1024 : : std::size_t ndofmax,
1025 : : tk::real tolref,
1026 : : std::vector< std::size_t >& ndofel ) const
1027 : : {
1028 : 0 : const auto& esuel = fd.Esuel();
1029 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
1030 : :
1031 [ - - ]: 0 : if(indicator == inciter::ctr::PrefIndicatorType::SPECTRAL_DECAY)
1032 : 0 : spectral_decay(nmat, nunk, esuel, unk, ndof, ndofmax, tolref, ndofel);
1033 : : else
1034 [ - - ][ - - ]: 0 : Throw( "No such adaptive indicator type" );
[ - - ]
1035 : 0 : }
1036 : :
1037 : : //! Compute the minimum time step size
1038 : : //! \param[in] fd Face connectivity and boundary conditions object
1039 : : //! \param[in] geoFace Face geometry array
1040 : : //! \param[in] geoElem Element geometry array
1041 : : // //! \param[in] ndofel Vector of local number of degrees of freedom
1042 : : //! \param[in] U Solution vector at recent time step
1043 : : //! \param[in] P Vector of primitive quantities at recent time step
1044 : : //! \param[in] nielem Number of internal elements
1045 : : //! \return Minimum time step size
1046 : : //! \details The allowable dt is calculated by looking at the maximum
1047 : : //! wave-speed in elements surrounding each face, times the area of that
1048 : : //! face. Once the maximum of this quantity over the mesh is determined,
1049 : : //! the volume of each cell is divided by this quantity. A minimum of this
1050 : : //! ratio is found over the entire mesh, which gives the allowable dt.
1051 : 260 : tk::real dt( const std::array< std::vector< tk::real >, 3 >&,
1052 : : const std::vector< std::size_t >&,
1053 : : const inciter::FaceData& fd,
1054 : : const tk::Fields& geoFace,
1055 : : const tk::Fields& geoElem,
1056 : : const std::vector< std::size_t >& /*ndofel*/,
1057 : : const tk::Fields& U,
1058 : : const tk::Fields& P,
1059 : : const std::size_t nielem ) const
1060 : : {
1061 : 260 : const auto ndof = g_inputdeck.get< tag::ndof >();
1062 : 260 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
1063 : :
1064 : 260 : auto mindt = timeStepSizeMultiMat( m_mat_blk, fd.Esuf(), geoFace, geoElem,
1065 : : nielem, nmat, U, P);
1066 : :
1067 : 260 : tk::real dgp = 0.0;
1068 [ + - ]: 260 : if (ndof == 4)
1069 : : {
1070 : 260 : dgp = 1.0;
1071 : : }
1072 [ - - ]: 0 : else if (ndof == 10)
1073 : : {
1074 : 0 : dgp = 2.0;
1075 : : }
1076 : :
1077 : : // Scale smallest dt with CFL coefficient and the CFL is scaled by (2*p+1)
1078 : : // where p is the order of the DG polynomial by linear stability theory.
1079 : 260 : mindt /= (2.0*dgp + 1.0);
1080 : 260 : return mindt;
1081 : : }
1082 : :
1083 : : //! Compute stiff terms for a single element
1084 : : //! \param[in] e Element number
1085 : : //! \param[in] geoElem Element geometry array
1086 : : //! \param[in] inpoel Element-node connectivity
1087 : : //! \param[in] coord Array of nodal coordinates
1088 : : //! \param[in] U Solution vector at recent time step
1089 : : //! \param[in] P Primitive vector at recent time step
1090 : : //! \param[in] ndofel Vector of local number of degrees of freedom
1091 : : //! \param[in,out] R Right-hand side vector computed
1092 : 0 : void stiff_rhs( std::size_t e,
1093 : : const tk::Fields& geoElem,
1094 : : const std::vector< std::size_t >& inpoel,
1095 : : const tk::UnsMesh::Coords& coord,
1096 : : const tk::Fields& U,
1097 : : const tk::Fields& P,
1098 : : const std::vector< std::size_t >& ndofel,
1099 : : tk::Fields& R ) const
1100 : : {
1101 : 0 : const auto ndof = g_inputdeck.get< tag::ndof >();
1102 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
1103 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
1104 : 0 : const auto intsharp =
1105 : 0 : g_inputdeck.get< tag::multimat, tag::intsharp >();
1106 : : const auto& solidx = inciter::g_inputdeck.get<
1107 : 0 : tag::matidxmap, tag::solidx >();
1108 : :
1109 [ - - ][ - - ]: 0 : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
1110 : : "vector and primitive vector at recent time step incorrect" );
1111 [ - - ][ - - ]: 0 : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
1112 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
1113 [ - - ][ - - ]: 0 : Assert( P.nprop() == rdof*m_nprim, "Number of components in primitive "
[ - - ][ - - ]
[ - - ]
1114 : : "vector must equal "+ std::to_string(rdof*m_nprim) );
1115 [ - - ][ - - ]: 0 : Assert( R.nprop() == ndof*nstiffeq(), "Number of components in "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
1116 : : "right-hand side must equal "+ std::to_string(ndof*nstiffeq()) );
1117 : :
1118 : : // set rhs to zero for element e
1119 [ - - ][ - - ]: 0 : for (std::size_t i=0; i<ndof*nstiffeq(); ++i)
1120 [ - - ]: 0 : R(e, i) = 0.0;
1121 : :
1122 : 0 : const auto& cx = coord[0];
1123 : 0 : const auto& cy = coord[1];
1124 : 0 : const auto& cz = coord[2];
1125 : :
1126 : 0 : auto ncomp = U.nprop()/rdof;
1127 : 0 : auto nprim = P.nprop()/rdof;
1128 : :
1129 [ - - ]: 0 : auto ng = tk::NGvol(ndofel[e]);
1130 : :
1131 : : // arrays for quadrature points
1132 : 0 : std::array< std::vector< tk::real >, 3 > coordgp;
1133 : 0 : std::vector< tk::real > wgp;
1134 : :
1135 [ - - ]: 0 : coordgp[0].resize( ng );
1136 [ - - ]: 0 : coordgp[1].resize( ng );
1137 [ - - ]: 0 : coordgp[2].resize( ng );
1138 [ - - ]: 0 : wgp.resize( ng );
1139 : :
1140 [ - - ]: 0 : tk::GaussQuadratureTet( ng, coordgp, wgp );
1141 : :
1142 : : // Extract the element coordinates
1143 : 0 : std::array< std::array< tk::real, 3>, 4 > coordel {{
1144 : : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
1145 : 0 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
1146 : 0 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
1147 : 0 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
1148 : : }};
1149 : :
1150 : : // Gaussian quadrature
1151 [ - - ]: 0 : for (std::size_t igp=0; igp<ng; ++igp)
1152 : : {
1153 : : // Compute the coordinates of quadrature point at physical domain
1154 [ - - ]: 0 : auto gp = tk::eval_gp( igp, coordel, coordgp );
1155 : :
1156 : : // Compute the basis function
1157 [ - - ]: 0 : auto B = tk::eval_basis( ndofel[e], coordgp[0][igp], coordgp[1][igp],
1158 : : coordgp[2][igp] );
1159 : :
1160 [ - - ]: 0 : auto state = tk::evalPolynomialSol(m_mat_blk, intsharp, ncomp, nprim,
1161 : : rdof, nmat, e, ndofel[e], inpoel, coord, geoElem, gp, B, U, P);
1162 : :
1163 : : // compute source
1164 : : // Loop through materials
1165 : 0 : std::size_t ksld = 0;
1166 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
1167 : : {
1168 [ - - ]: 0 : if (solidx[k] > 0)
1169 : : {
1170 : 0 : tk::real alpha = state[inciter::volfracIdx(nmat, k)];
1171 : : std::array< std::array< tk::real, 3 >, 3 > g;
1172 : : // Compute the source terms
1173 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1174 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
1175 : 0 : g[i][j] = state[inciter::deformIdx(nmat,solidx[k],i,j)];
1176 : :
1177 : : // Compute Lp
1178 : : // Reference: Ortega, A. L., Lombardini, M., Pullin, D. I., &
1179 : : // Meiron, D. I. (2014). Numerical simulation of elastic–plastic
1180 : : // solid mechanics using an Eulerian stretch tensor approach and
1181 : : // HLLD Riemann solver. Journal of Computational Physics, 257,
1182 : : // 414-441
1183 : : std::array< std::array< tk::real, 3 >, 3 > Lp;
1184 : :
1185 : : // 1. Compute dev(sigma)
1186 : 0 : auto sigma_dev = m_mat_blk[k].computeTensor< EOS::CauchyStress >(
1187 [ - - ]: 0 : 0.0, 0.0, 0.0, 0.0, 0.0, alpha, k, g );
1188 : 0 : tk::real apr = state[ncomp+inciter::pressureIdx(nmat, k)];
1189 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) sigma_dev[i][i] -= apr;
1190 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1191 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
1192 : 0 : sigma_dev[i][j] /= alpha;
1193 : 0 : tk::real sigma_trace =
1194 : 0 : sigma_dev[0][0]+sigma_dev[1][1]+sigma_dev[2][2];
1195 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1196 : 0 : sigma_dev[i][i] -= sigma_trace/3.0;
1197 : :
1198 : : // 2. Compute inv(g)
1199 : : double ginv[9];
1200 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1201 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
1202 : 0 : ginv[3*i+j] = g[i][j];
1203 : : lapack_int ipiv[3];
1204 : : #ifndef NDEBUG
1205 : : lapack_int ierr =
1206 : : #endif
1207 [ - - ]: 0 : LAPACKE_dgetrf(LAPACK_ROW_MAJOR, 3, 3, ginv, 3, ipiv);
1208 [ - - ][ - - ]: 0 : Assert(ierr==0, "Lapack error in LU factorization of g");
[ - - ][ - - ]
1209 : : #ifndef NDEBUG
1210 : : lapack_int jerr =
1211 : : #endif
1212 [ - - ]: 0 : LAPACKE_dgetri(LAPACK_ROW_MAJOR, 3, ginv, 3, ipiv);
1213 [ - - ][ - - ]: 0 : Assert(jerr==0, "Lapack error in inverting g");
[ - - ][ - - ]
1214 : :
1215 : : // 3. Compute dev(sigma)*inv(g)
1216 : : std::array< std::array< tk::real, 3 >, 3 > aux_mat;
1217 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1218 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
1219 : : {
1220 : 0 : tk::real sum = 0.0;
1221 [ - - ]: 0 : for (std::size_t l=0; l<3; ++l)
1222 : 0 : sum += sigma_dev[i][l]*ginv[3*l+j];
1223 : 0 : aux_mat[i][j] = sum;
1224 : : }
1225 : :
1226 : : // 4. Compute g*(dev(sigma)*inv(g))
1227 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1228 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
1229 : : {
1230 : 0 : tk::real sum = 0.0;
1231 [ - - ]: 0 : for (std::size_t l=0; l<3; ++l)
1232 : 0 : sum += g[i][l]*aux_mat[l][j];
1233 : 0 : Lp[i][j] = sum;
1234 : : }
1235 : :
1236 : : // 5. Divide by 2*mu*tau
1237 : : // 'Perfect' plasticity
1238 [ - - ]: 0 : tk::real yield_stress = getmatprop< tag::yield_stress >(k);
1239 : 0 : tk::real equiv_stress = 0.0;
1240 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1241 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
1242 : 0 : equiv_stress += sigma_dev[i][j]*sigma_dev[i][j];
1243 : 0 : equiv_stress = std::sqrt(3.0*equiv_stress/2.0);
1244 : : // rel_factor = 1/tau <- Perfect plasticity for now.
1245 : 0 : tk::real rel_factor = 0.0;
1246 [ - - ]: 0 : if (equiv_stress >= yield_stress)
1247 : 0 : rel_factor = 1.0e07;
1248 [ - - ]: 0 : tk::real mu = getmatprop< tag::mu >(k);
1249 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1250 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
1251 : 0 : Lp[i][j] *= rel_factor/(2.0*mu);
1252 : :
1253 : : // Compute the source terms
1254 [ - - ]: 0 : std::vector< tk::real > s(9*ndof, 0.0);
1255 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1256 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
1257 [ - - ]: 0 : for (std::size_t idof=0; idof<ndof; ++idof)
1258 : : {
1259 : 0 : s[(i*3+j)*ndof+idof] = B[idof] * (Lp[i][0]*g[0][j]
1260 : 0 : +Lp[i][1]*g[1][j]
1261 : 0 : +Lp[i][2]*g[2][j]);
1262 : : }
1263 : :
1264 [ - - ]: 0 : auto wt = wgp[igp] * geoElem(e, 0);
1265 : :
1266 : : // Contribute to the right-hand-side
1267 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
1268 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
1269 [ - - ]: 0 : for (std::size_t idof=0; idof<ndof; ++idof)
1270 : : {
1271 : 0 : std::size_t srcId = (i*3+j)*ndof+idof;
1272 : 0 : std::size_t dofId = solidTensorIdx(ksld,i,j)*ndof+idof;
1273 [ - - ]: 0 : R(e, dofId) += wt * s[srcId];
1274 : : }
1275 : :
1276 : 0 : ksld++;
1277 : : }
1278 : : }
1279 : :
1280 : : }
1281 : 0 : }
1282 : :
1283 : : //! Extract the velocity field at cell nodes. Currently unused.
1284 : : //! \param[in] U Solution vector at recent time step
1285 : : //! \param[in] N Element node indices
1286 : : //! \return Array of the four values of the velocity field
1287 : : std::array< std::array< tk::real, 4 >, 3 >
1288 : : velocity( const tk::Fields& U,
1289 : : const std::array< std::vector< tk::real >, 3 >&,
1290 : : const std::array< std::size_t, 4 >& N ) const
1291 : : {
1292 : : const auto rdof = g_inputdeck.get< tag::rdof >();
1293 : : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
1294 : :
1295 : : std::array< std::array< tk::real, 4 >, 3 > v;
1296 : : v[0] = U.extract( momentumDofIdx(nmat, 0, rdof, 0), N );
1297 : : v[1] = U.extract( momentumDofIdx(nmat, 1, rdof, 0), N );
1298 : : v[2] = U.extract( momentumDofIdx(nmat, 2, rdof, 0), N );
1299 : :
1300 : : std::vector< std::array< tk::real, 4 > > ar;
1301 : : ar.resize(nmat);
1302 : : for (std::size_t k=0; k<nmat; ++k)
1303 : : ar[k] = U.extract( densityDofIdx(nmat, k, rdof, 0), N );
1304 : :
1305 : : std::array< tk::real, 4 > r{{ 0.0, 0.0, 0.0, 0.0 }};
1306 : : for (std::size_t i=0; i<r.size(); ++i) {
1307 : : for (std::size_t k=0; k<nmat; ++k)
1308 : : r[i] += ar[k][i];
1309 : : }
1310 : :
1311 : : std::transform( r.begin(), r.end(), v[0].begin(), v[0].begin(),
1312 : : []( tk::real s, tk::real& d ){ return d /= s; } );
1313 : : std::transform( r.begin(), r.end(), v[1].begin(), v[1].begin(),
1314 : : []( tk::real s, tk::real& d ){ return d /= s; } );
1315 : : std::transform( r.begin(), r.end(), v[2].begin(), v[2].begin(),
1316 : : []( tk::real s, tk::real& d ){ return d /= s; } );
1317 : : return v;
1318 : : }
1319 : :
1320 : : //! Return a map that associates user-specified strings to functions
1321 : : //! \return Map that associates user-specified strings to functions that
1322 : : //! compute relevant quantities to be output to file
1323 : 310 : std::map< std::string, tk::GetVarFn > OutVarFn() const
1324 : 310 : { return MultiMatOutVarFn(); }
1325 : :
1326 : : //! Return analytic field names to be output to file
1327 : : //! \return Vector of strings labelling analytic fields output in file
1328 : 0 : std::vector< std::string > analyticFieldNames() const {
1329 : 0 : auto nmat = g_inputdeck.get< eq, tag::nmat >();
1330 : :
1331 : 0 : return MultiMatFieldNames(nmat);
1332 : : }
1333 : :
1334 : : //! Return time history field names to be output to file
1335 : : //! \return Vector of strings labelling time history fields output in file
1336 : 0 : std::vector< std::string > histNames() const {
1337 : 0 : return MultiMatHistNames();
1338 : : }
1339 : :
1340 : : //! Return surface field output going to file
1341 : : std::vector< std::vector< tk::real > >
1342 : 0 : surfOutput( const std::map< int, std::vector< std::size_t > >&,
1343 : : tk::Fields& ) const
1344 : : {
1345 : 0 : std::vector< std::vector< tk::real > > s; // punt for now
1346 : 0 : return s;
1347 : : }
1348 : :
1349 : : //! Return time history field output evaluated at time history points
1350 : : //! \param[in] h History point data
1351 : : //! \param[in] inpoel Element-node connectivity
1352 : : //! \param[in] coord Array of nodal coordinates
1353 : : //! \param[in] U Array of unknowns
1354 : : //! \param[in] P Array of primitive quantities
1355 : : //! \return Vector of time history output of bulk flow quantities (density,
1356 : : //! velocity, total energy, and pressure) evaluated at time history points
1357 : : std::vector< std::vector< tk::real > >
1358 : 0 : histOutput( const std::vector< HistData >& h,
1359 : : const std::vector< std::size_t >& inpoel,
1360 : : const tk::UnsMesh::Coords& coord,
1361 : : const tk::Fields& U,
1362 : : const tk::Fields& P ) const
1363 : : {
1364 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
1365 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
1366 : :
1367 : 0 : const auto& x = coord[0];
1368 : 0 : const auto& y = coord[1];
1369 : 0 : const auto& z = coord[2];
1370 : :
1371 [ - - ]: 0 : std::vector< std::vector< tk::real > > Up(h.size());
1372 : :
1373 : 0 : std::size_t j = 0;
1374 [ - - ]: 0 : for (const auto& p : h) {
1375 : 0 : auto e = p.get< tag::elem >();
1376 : 0 : auto chp = p.get< tag::coord >();
1377 : :
1378 : : // Evaluate inverse Jacobian
1379 : 0 : std::array< std::array< tk::real, 3>, 4 > cp{{
1380 : : {{ x[inpoel[4*e ]], y[inpoel[4*e ]], z[inpoel[4*e ]] }},
1381 : 0 : {{ x[inpoel[4*e+1]], y[inpoel[4*e+1]], z[inpoel[4*e+1]] }},
1382 : 0 : {{ x[inpoel[4*e+2]], y[inpoel[4*e+2]], z[inpoel[4*e+2]] }},
1383 : 0 : {{ x[inpoel[4*e+3]], y[inpoel[4*e+3]], z[inpoel[4*e+3]] }} }};
1384 : 0 : auto J = tk::inverseJacobian( cp[0], cp[1], cp[2], cp[3] );
1385 : :
1386 : : // evaluate solution at history-point
1387 : 0 : std::array< tk::real, 3 > dc{{chp[0]-cp[0][0], chp[1]-cp[0][1],
1388 : 0 : chp[2]-cp[0][2]}};
1389 [ - - ]: 0 : auto B = tk::eval_basis(rdof, tk::dot(J[0],dc), tk::dot(J[1],dc),
1390 : 0 : tk::dot(J[2],dc));
1391 [ - - ]: 0 : auto uhp = eval_state(m_ncomp, rdof, rdof, e, U, B);
1392 [ - - ]: 0 : auto php = eval_state(m_nprim, rdof, rdof, e, P, B);
1393 : :
1394 : : // store solution in history output vector
1395 [ - - ]: 0 : Up[j].resize(6, 0.0);
1396 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
1397 : 0 : Up[j][0] += uhp[densityIdx(nmat,k)];
1398 : 0 : Up[j][4] += uhp[energyIdx(nmat,k)];
1399 : 0 : Up[j][5] += php[pressureIdx(nmat,k)];
1400 : : }
1401 : 0 : Up[j][1] = php[velocityIdx(nmat,0)];
1402 : 0 : Up[j][2] = php[velocityIdx(nmat,1)];
1403 : 0 : Up[j][3] = php[velocityIdx(nmat,2)];
1404 : 0 : ++j;
1405 : : }
1406 : :
1407 : 0 : return Up;
1408 : : }
1409 : :
1410 : : //! Return names of integral variables to be output to diagnostics file
1411 : : //! \return Vector of strings labelling integral variables output
1412 : 12 : std::vector< std::string > names() const
1413 : : {
1414 : 12 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
1415 : 12 : return MultiMatDiagNames(nmat);
1416 : : }
1417 : :
1418 : : //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
1419 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
1420 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
1421 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
1422 : : //! \param[in] t Physical time at which to evaluate the analytic solution
1423 : : //! \return Vector of analytic solution at given location and time
1424 : : std::vector< tk::real >
1425 : 0 : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
1426 : 0 : { return Problem::analyticSolution( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
1427 : :
1428 : : //! Return analytic solution for conserved variables
1429 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
1430 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
1431 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
1432 : : //! \param[in] t Physical time at which to evaluate the analytic solution
1433 : : //! \return Vector of analytic solution at given location and time
1434 : : std::vector< tk::real >
1435 : 1322920 : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
1436 : 1322920 : { return Problem::initialize( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
1437 : :
1438 : : //! Return cell-averaged specific total energy for an element
1439 : : //! \param[in] e Element id for which total energy is required
1440 : : //! \param[in] unk Vector of conserved quantities
1441 : : //! \return Cell-averaged specific total energy for given element
1442 : 822640 : tk::real sp_totalenergy(std::size_t e, const tk::Fields& unk) const
1443 : : {
1444 : 822640 : const auto rdof = g_inputdeck.get< tag::rdof >();
1445 : 822640 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
1446 : :
1447 : 822640 : tk::real sp_te(0.0);
1448 : : // sum each material total energy
1449 [ + + ]: 2858500 : for (std::size_t k=0; k<nmat; ++k) {
1450 : 2035860 : sp_te += unk(e, energyDofIdx(nmat,k,rdof,0));
1451 : : }
1452 : 822640 : return sp_te;
1453 : : }
1454 : :
1455 : : private:
1456 : : //! Number of components in this PDE system
1457 : : const ncomp_t m_ncomp;
1458 : : //! Number of primitive quantities stored in this PDE system
1459 : : const ncomp_t m_nprim;
1460 : : //! Riemann solver
1461 : : tk::RiemannFluxFn m_riemann;
1462 : : //! BC configuration
1463 : : BCStateFn m_bc;
1464 : : //! EOS material block
1465 : : std::vector< EOS > m_mat_blk;
1466 : :
1467 : : //! Evaluate conservative part of physical flux function for this PDE system
1468 : : //! \param[in] ncomp Number of scalar components in this PDE system
1469 : : //! \param[in] ugp Numerical solution at the Gauss point at which to
1470 : : //! evaluate the flux
1471 : : //! \return Flux vectors for all components in this PDE system
1472 : : //! \note The function signature must follow tk::FluxFn
1473 : : static tk::FluxFn::result_type
1474 : 7374000 : flux( ncomp_t ncomp,
1475 : : const std::vector< EOS >& mat_blk,
1476 : : const std::vector< tk::real >& ugp,
1477 : : const std::vector< std::array< tk::real, 3 > >& )
1478 : : {
1479 : 7374000 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
1480 : :
1481 : 7374000 : return tk::fluxTerms(ncomp, nmat, mat_blk, ugp);
1482 : : }
1483 : :
1484 : : //! \brief Boundary state function providing the left and right state of a
1485 : : //! face at Dirichlet boundaries
1486 : : //! \param[in] ncomp Number of scalar components in this PDE system
1487 : : //! \param[in] mat_blk EOS material block
1488 : : //! \param[in] ul Left (domain-internal) state
1489 : : //! \param[in] x X-coordinate at which to compute the states
1490 : : //! \param[in] y Y-coordinate at which to compute the states
1491 : : //! \param[in] z Z-coordinate at which to compute the states
1492 : : //! \param[in] t Physical time
1493 : : //! \return Left and right states for all scalar components in this PDE
1494 : : //! system
1495 : : //! \note The function signature must follow tk::StateFn. For multimat, the
1496 : : //! left or right state is the vector of conserved quantities, followed by
1497 : : //! the vector of primitive quantities appended to it.
1498 : : static tk::StateFn::result_type
1499 : 48000 : dirichlet( ncomp_t ncomp,
1500 : : const std::vector< EOS >& mat_blk,
1501 : : const std::vector< tk::real >& ul, tk::real x, tk::real y,
1502 : : tk::real z, tk::real t, const std::array< tk::real, 3 >& )
1503 : : {
1504 : 48000 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
1505 : : const auto& solidx = g_inputdeck.get<
1506 : 48000 : tag::matidxmap, tag::solidx >();
1507 : :
1508 [ + - ]: 48000 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
1509 : :
1510 [ + - ]: 96000 : auto ur = Problem::initialize( ncomp, mat_blk, x, y, z, t );
1511 [ - + ][ - - ]: 48000 : Assert( ur.size() == ncomp, "Incorrect size for boundary state vector" );
[ - - ][ - - ]
1512 : :
1513 [ + - ]: 48000 : ur.resize(ul.size());
1514 : :
1515 : 48000 : tk::real rho(0.0);
1516 [ + + ]: 192000 : for (std::size_t k=0; k<nmat; ++k)
1517 : 144000 : rho += ur[densityIdx(nmat, k)];
1518 : :
1519 : : // get primitives in boundary state
1520 : :
1521 : : // velocity
1522 : 48000 : ur[ncomp+velocityIdx(nmat, 0)] = ur[momentumIdx(nmat, 0)] / rho;
1523 : 48000 : ur[ncomp+velocityIdx(nmat, 1)] = ur[momentumIdx(nmat, 1)] / rho;
1524 : 48000 : ur[ncomp+velocityIdx(nmat, 2)] = ur[momentumIdx(nmat, 2)] / rho;
1525 : :
1526 : : // material pressures
1527 [ + + ]: 192000 : for (std::size_t k=0; k<nmat; ++k)
1528 : : {
1529 [ + - ]: 144000 : auto gk = getDeformGrad(nmat, k, ur);
1530 [ + - ]: 288000 : ur[ncomp+pressureIdx(nmat, k)] = mat_blk[k].compute< EOS::pressure >(
1531 : 144000 : ur[densityIdx(nmat, k)], ur[ncomp+velocityIdx(nmat, 0)],
1532 : 144000 : ur[ncomp+velocityIdx(nmat, 1)], ur[ncomp+velocityIdx(nmat, 2)],
1533 : 144000 : ur[energyIdx(nmat, k)], ur[volfracIdx(nmat, k)], k, gk );
1534 : : }
1535 : :
1536 [ - + ][ - - ]: 48000 : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
1537 : : "boundary state vector" );
1538 : :
1539 [ + - ]: 96000 : return {{ std::move(ul), std::move(ur) }};
1540 : : }
1541 : :
1542 : : // Other boundary condition types that do not depend on "Problem" should be
1543 : : // added in BCFunctions.hpp
1544 : : };
1545 : :
1546 : : } // dg::
1547 : :
1548 : : } // inciter::
1549 : :
1550 : : #endif // DGMultiMat_h
|
