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