Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/MultiMat/FVMultiMat.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 finite volumes
9 : : \details This file implements calls to the physics operators governing
10 : : compressible multi-material flow (with velocity equilibrium) using finite
11 : : volume discretizations.
12 : : */
13 : : // *****************************************************************************
14 : : #ifndef FVMultiMat_h
15 : : #define FVMultiMat_h
16 : :
17 : : #include <cmath>
18 : : #include <algorithm>
19 : : #include <unordered_set>
20 : : #include <map>
21 : : #include <array>
22 : :
23 : : #include "Macro.hpp"
24 : : #include "Exception.hpp"
25 : : #include "Vector.hpp"
26 : : #include "ContainerUtil.hpp"
27 : : #include "UnsMesh.hpp"
28 : : #include "Inciter/InputDeck/InputDeck.hpp"
29 : : #include "Integrate/Basis.hpp"
30 : : #include "Integrate/Quadrature.hpp"
31 : : #include "Integrate/Initialize.hpp"
32 : : #include "Integrate/Mass.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 "RiemannChoice.hpp"
39 : : #include "MultiMat/MultiMatIndexing.hpp"
40 : : #include "Reconstruction.hpp"
41 : : #include "Limiter.hpp"
42 : : #include "Problem/FieldOutput.hpp"
43 : : #include "Problem/BoxInitialization.hpp"
44 : : #include "MultiMat/BCFunctions.hpp"
45 : : #include "MultiMat/MiscMultiMatFns.hpp"
46 : :
47 : : namespace inciter {
48 : :
49 : : extern ctr::InputDeck g_inputdeck;
50 : :
51 : : namespace fv {
52 : :
53 : : //! \brief MultiMat used polymorphically with tk::FVPDE
54 : : //! \details The template arguments specify policies and are used to configure
55 : : //! the behavior of the class. The policies are:
56 : : //! - Physics - physics configuration, see PDE/MultiMat/Physics.h
57 : : //! - Problem - problem configuration, see PDE/MultiMat/Problem.h
58 : : //! \note The default physics is Euler, set in inciter::deck::check_multimat()
59 : : template< class Physics, class Problem >
60 : : class MultiMat {
61 : :
62 : : private:
63 : : using eq = tag::multimat;
64 : :
65 : : public:
66 : : //! Constructor
67 : 79 : explicit MultiMat() :
68 : : m_physics(),
69 : 79 : m_ncomp( g_inputdeck.get< tag::ncomp >() ),
70 : : m_riemann( multimatRiemannSolver(
71 : 79 : g_inputdeck.get< tag::flux >() ) )
72 : : {
73 : : // associate boundary condition configurations with state functions
74 [ + - ][ + - ]: 553 : brigand::for_each< ctr::bclist::Keys >( ConfigBC( m_bc,
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ][ - - ]
75 : : { dirichlet
76 : : , symmetry
77 : : , invalidBC // Inlet BC not implemented
78 : : , invalidBC // Outlet BC not implemented
79 : : , farfieldOutlet
80 : : , extrapolate } ) );
81 : :
82 : : // EoS initialization
83 [ + - ]: 79 : initializeMaterialEoS( m_mat_blk );
84 : 79 : }
85 : :
86 : : //! Find the number of primitive quantities required for this PDE system
87 : : //! \return The number of primitive quantities required to be stored for
88 : : //! this PDE system
89 : 16324 : std::size_t nprim() const
90 : : {
91 : 16324 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
92 : : // multimat needs individual material pressures and velocities currently
93 : 16324 : return (nmat+3);
94 : : }
95 : :
96 : : //! Find the number of materials set up for this PDE system
97 : : //! \return The number of materials set up for this PDE system
98 : 0 : std::size_t nmat() const
99 : : {
100 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
101 : 0 : return nmat;
102 : : }
103 : :
104 : : //! Determine elements that lie inside the user-defined IC box
105 : : //! \param[in] geoElem Element geometry array
106 : : //! \param[in] nielem Number of internal elements
107 : : //! \param[in,out] inbox List of nodes at which box user ICs are set for
108 : : //! each IC box
109 : 297 : void IcBoxElems( const tk::Fields& geoElem,
110 : : std::size_t nielem,
111 : : std::vector< std::unordered_set< std::size_t > >& inbox ) const
112 : : {
113 : 297 : tk::BoxElems< eq >(geoElem, nielem, inbox);
114 : 297 : }
115 : :
116 : : //! Initalize the compressible flow equations, prepare for time integration
117 : : //! \param[in] L Block diagonal mass matrix
118 : : //! \param[in] inpoel Element-node connectivity
119 : : //! \param[in] coord Array of nodal coordinates
120 : : //! \param[in] inbox List of elements at which box user ICs are set for
121 : : //! each IC box
122 : : //! \param[in] elemblkid Element ids associated with mesh block ids where
123 : : //! user ICs are set
124 : : //! \param[in,out] unk Array of unknowns
125 : : //! \param[in] t Physical time
126 : : //! \param[in] nielem Number of internal elements
127 : 297 : void initialize( const tk::Fields& L,
128 : : const std::vector< std::size_t >& inpoel,
129 : : const tk::UnsMesh::Coords& coord,
130 : : const std::vector< std::unordered_set< std::size_t > >& inbox,
131 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
132 : : elemblkid,
133 : : tk::Fields& unk,
134 : : tk::real t,
135 : : const std::size_t nielem ) const
136 : : {
137 [ + - ][ + - ]: 297 : tk::initialize( m_ncomp, m_mat_blk, L, inpoel, coord,
138 : : Problem::initialize, unk, t, nielem );
139 : :
140 : 297 : const auto rdof = g_inputdeck.get< tag::rdof >();
141 : 297 : const auto& ic = g_inputdeck.get< tag::ic >();
142 : 297 : const auto& icbox = ic.get< tag::box >();
143 : 297 : const auto& icmbk = ic.get< tag::meshblock >();
144 : :
145 : 297 : const auto& bgpre = ic.get< tag::pressure >();
146 : 297 : const auto& bgtemp = ic.get< tag::temperature >();
147 : :
148 : : // Set initial conditions inside user-defined IC boxes and mesh blocks
149 [ + - ]: 594 : std::vector< tk::real > s(m_ncomp, 0.0);
150 [ + + ]: 15877 : for (std::size_t e=0; e<nielem; ++e) {
151 : : // inside user-defined box
152 [ - + ]: 15580 : if (!icbox.empty()) {
153 : 0 : std::size_t bcnt = 0;
154 [ - - ]: 0 : for (const auto& b : icbox) { // for all boxes
155 [ - - ][ - - ]: 0 : if (inbox.size() > bcnt && inbox[bcnt].find(e) != inbox[bcnt].end())
[ - - ][ - - ]
156 : : {
157 [ - - ]: 0 : std::vector< tk::real > box
158 : 0 : { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
159 : 0 : b.template get< tag::ymin >(), b.template get< tag::ymax >(),
160 : 0 : b.template get< tag::zmin >(), b.template get< tag::zmax >() };
161 : 0 : auto V_ex = (box[1]-box[0]) * (box[3]-box[2]) * (box[5]-box[4]);
162 [ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) {
163 : 0 : auto mark = c*rdof;
164 [ - - ]: 0 : s[c] = unk(e,mark);
165 : : // set high-order DOFs to zero
166 [ - - ]: 0 : for (std::size_t i=1; i<rdof; ++i)
167 [ - - ]: 0 : unk(e,mark+i) = 0.0;
168 : : }
169 [ - - ]: 0 : initializeBox<ctr::boxList>( m_mat_blk, V_ex, t, b, bgpre,
170 : : bgtemp, s );
171 : : // store box-initialization in solution vector
172 [ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) {
173 : 0 : auto mark = c*rdof;
174 [ - - ]: 0 : unk(e,mark) = s[c];
175 : : }
176 : : }
177 : 0 : ++bcnt;
178 : : }
179 : : }
180 : :
181 : : // inside user-specified mesh blocks
182 [ + + ]: 16504 : for (const auto& b : icmbk) { // for all blocks
183 : 924 : auto blid = b.get< tag::blockid >();
184 : 924 : auto V_ex = b.get< tag::volume >();
185 [ + - ][ + - ]: 924 : if (elemblkid.find(blid) != elemblkid.end()) {
186 [ + - ]: 924 : const auto& elset = tk::cref_find(elemblkid, blid);
187 [ + - ][ + + ]: 924 : if (elset.find(e) != elset.end()) {
188 [ + - ]: 104 : initializeBox<ctr::meshblockList>( m_mat_blk, V_ex, t, b,
189 : : bgpre, bgtemp, s );
190 : : // store initialization in solution vector
191 [ + + ]: 1040 : for (std::size_t c=0; c<m_ncomp; ++c) {
192 : 936 : auto mark = c*rdof;
193 [ + - ]: 936 : unk(e,mark) = s[c];
194 : : }
195 : : }
196 : : }
197 : : }
198 : : }
199 : 297 : }
200 : :
201 : : //! Compute the left hand side block-diagonal mass matrix
202 : : //! \param[in] geoElem Element geometry array
203 : : //! \param[in,out] l Block diagonal mass matrix
204 : 297 : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const {
205 : 297 : const auto nelem = geoElem.nunk();
206 [ + + ]: 53576 : for (std::size_t e=0; e<nelem; ++e)
207 [ + + ]: 679148 : for (ncomp_t c=0; c<m_ncomp; ++c)
208 : 625869 : l(e, c) = geoElem(e,0);
209 : 297 : }
210 : :
211 : : //! Update the primitives for this PDE system
212 : : //! \param[in] unk Array of unknowns
213 : : //! \param[in,out] prim Array of primitives
214 : : //! \param[in] nielem Number of internal elements
215 : : //! \details This function computes and stores the dofs for primitive
216 : : //! quantities, which are required for obtaining reconstructed states used
217 : : //! in the Riemann solver. See /PDE/Riemann/AUSM.hpp, where the
218 : : //! normal velocity for advection is calculated from independently
219 : : //! reconstructed velocities.
220 : 1915 : void updatePrimitives( const tk::Fields& unk,
221 : : tk::Fields& prim,
222 : : std::size_t nielem ) const
223 : : {
224 : 1915 : const auto rdof = g_inputdeck.get< tag::rdof >();
225 : 1915 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
226 : :
227 [ - + ][ - - ]: 1915 : Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
228 : : "vector and primitive vector at recent time step incorrect" );
229 [ - + ][ - - ]: 1915 : Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
230 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
231 [ - + ][ - - ]: 1915 : Assert( prim.nprop() == rdof*nprim(), "Number of components in vector of "
[ - - ][ - - ]
[ - - ]
232 : : "primitive quantities must equal "+ std::to_string(rdof*nprim()) );
233 : :
234 [ + + ]: 238255 : for (std::size_t e=0; e<nielem; ++e)
235 : : {
236 : : // cell-average bulk density
237 : 236340 : tk::real rhob(0.0);
238 [ + + ]: 774720 : for (std::size_t k=0; k<nmat; ++k)
239 : : {
240 [ + - ]: 538380 : rhob += unk(e, densityDofIdx(nmat, k, rdof, 0));
241 : : }
242 : :
243 : : // cell-average velocity
244 : : std::array< tk::real, 3 >
245 [ + - ]: 236340 : vel{{ unk(e, momentumDofIdx(nmat, 0, rdof, 0))/rhob,
246 [ + - ]: 236340 : unk(e, momentumDofIdx(nmat, 1, rdof, 0))/rhob,
247 [ + - ]: 236340 : unk(e, momentumDofIdx(nmat, 2, rdof, 0))/rhob }};
248 : :
249 [ + + ]: 945360 : for (std::size_t idir=0; idir<3; ++idir)
250 : : {
251 [ + - ]: 709020 : prim(e, velocityDofIdx(nmat, idir, rdof, 0)) = vel[idir];
252 [ + + ]: 2836080 : for (std::size_t idof=1; idof<rdof; ++idof)
253 [ + - ]: 2127060 : prim(e, velocityDofIdx(nmat, idir, rdof, idof)) = 0.0;
254 : : }
255 : :
256 : : // cell-average material pressure
257 [ + + ]: 774720 : for (std::size_t k=0; k<nmat; ++k)
258 : : {
259 [ + - ]: 538380 : tk::real arhomat = unk(e, densityDofIdx(nmat, k, rdof, 0));
260 [ + - ]: 538380 : tk::real arhoemat = unk(e, energyDofIdx(nmat, k, rdof, 0));
261 [ + - ]: 538380 : tk::real alphamat = unk(e, volfracDofIdx(nmat, k, rdof, 0));
262 [ + - ][ + - ]: 538380 : auto gmat = getDeformGrad(nmat, k, unk.extract(e));
263 [ + - ]: 538380 : prim(e, pressureDofIdx(nmat, k, rdof, 0)) =
264 [ + - ]: 538380 : m_mat_blk[k].compute< EOS::pressure >( arhomat, vel[0], vel[1],
265 : : vel[2], arhoemat, alphamat, k, gmat );
266 [ + - ]: 538380 : prim(e, pressureDofIdx(nmat, k, rdof, 0)) =
267 : 538380 : constrain_pressure( m_mat_blk,
268 [ + - ][ + - ]: 538380 : prim(e, pressureDofIdx(nmat, k, rdof, 0)), arhomat, alphamat, k);
269 [ + + ]: 2153520 : for (std::size_t idof=1; idof<rdof; ++idof)
270 [ + - ]: 1615140 : prim(e, pressureDofIdx(nmat, k, rdof, idof)) = 0.0;
271 : : }
272 : : }
273 : 1915 : }
274 : :
275 : : //! Clean up the state of trace materials for this PDE system
276 : : //! \param[in] t Physical time
277 : : //! \param[in] geoElem Element geometry array
278 : : //! \param[in,out] unk Array of unknowns
279 : : //! \param[in,out] prim Array of primitives
280 : : //! \param[in] nielem Number of internal elements
281 : : //! \details This function cleans up the state of materials present in trace
282 : : //! quantities in each cell. Specifically, the state of materials with
283 : : //! very low volume-fractions in a cell is replaced by the state of the
284 : : //! material which is present in the largest quantity in that cell. This
285 : : //! becomes necessary when shocks pass through cells which contain a very
286 : : //! small amount of material. The state of that tiny material might
287 : : //! become unphysical and cause solution to diverge; thus requiring such
288 : : //! a "reset".
289 : 1618 : void cleanTraceMaterial( tk::real t,
290 : : const tk::Fields& geoElem,
291 : : tk::Fields& unk,
292 : : tk::Fields& prim,
293 : : std::size_t nielem ) const
294 : : {
295 : 1618 : [[maybe_unused]] const auto rdof = g_inputdeck.get< tag::rdof >();
296 : 1618 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
297 : :
298 [ - + ][ - - ]: 1618 : Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
299 : : "vector and primitive vector at recent time step incorrect" );
300 [ - + ][ - - ]: 1618 : Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
301 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
302 [ - + ][ - - ]: 1618 : Assert( prim.nprop() == rdof*nprim(), "Number of components in vector of "
[ - - ][ - - ]
[ - - ]
303 : : "primitive quantities must equal "+ std::to_string(rdof*nprim()) );
304 [ - + ][ - - ]: 1618 : Assert( (g_inputdeck.get< tag::ndof >()) <= 4, "High-order "
[ - - ][ - - ]
305 : : "discretizations not set up for multimat cleanTraceMaterial()" );
306 : :
307 : 1618 : auto neg_density = cleanTraceMultiMat(t, nielem, m_mat_blk, geoElem, nmat,
308 : : unk, prim);
309 : :
310 [ - + ][ - - ]: 1618 : if (neg_density) Throw("Negative partial density.");
[ - - ][ - - ]
311 : 1618 : }
312 : :
313 : : //! Reconstruct second-order solution from first-order
314 : : //! \param[in] geoElem Element geometry array
315 : : //! \param[in] fd Face connectivity and boundary conditions object
316 : : //! \param[in] esup Elements-surrounding-nodes connectivity
317 : : //! \param[in] inpoel Element-node connectivity
318 : : //! \param[in] coord Array of nodal coordinates
319 : : //! \param[in,out] U Solution vector at recent time step
320 : : //! \param[in,out] P Vector of primitives at recent time step
321 : 1618 : void reconstruct( const tk::Fields& geoElem,
322 : : const inciter::FaceData& fd,
323 : : const std::map< std::size_t, std::vector< std::size_t > >&
324 : : esup,
325 : : const std::vector< std::size_t >& inpoel,
326 : : const tk::UnsMesh::Coords& coord,
327 : : tk::Fields& U,
328 : : tk::Fields& P ) const
329 : : {
330 : 1618 : const auto rdof = g_inputdeck.get< tag::rdof >();
331 : 1618 : const auto nelem = fd.Esuel().size()/4;
332 : 1618 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
333 : :
334 [ - + ][ - - ]: 1618 : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
335 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
336 : :
337 : : //----- reconstruction of conserved quantities -----
338 : : //--------------------------------------------------
339 : : // specify how many variables need to be reconstructed
340 : 3236 : std::vector< std::size_t > vars;
341 [ + + ]: 6022 : for (std::size_t k=0; k<nmat; ++k) {
342 [ + - ]: 4404 : vars.push_back(volfracIdx(nmat,k));
343 [ + - ]: 4404 : vars.push_back(densityIdx(nmat,k));
344 : : }
345 : :
346 : : // 1. solve 3x3 least-squares system
347 [ + + ]: 222378 : for (std::size_t e=0; e<nelem; ++e)
348 : : {
349 : : // Reconstruct second-order dofs of volume-fractions in Taylor space
350 : : // using nodal-stencils, for a good interface-normal estimate
351 [ + - ]: 220760 : tk::recoLeastSqExtStencil( rdof, e, esup, inpoel, geoElem, U, vars );
352 : : }
353 : :
354 : : // 2. transform reconstructed derivatives to Dubiner dofs
355 [ + - ]: 1618 : tk::transform_P0P1(rdof, nelem, inpoel, coord, U, vars);
356 : :
357 : : //----- reconstruction of primitive quantities -----
358 : : //--------------------------------------------------
359 : : // For multimat, conserved and primitive quantities are reconstructed
360 : : // separately.
361 : 1618 : vars.clear();
362 [ + + ][ + - ]: 10876 : for (std::size_t c=0; c<nprim(); ++c) vars.push_back(c);
363 : : // 1.
364 [ + + ]: 222378 : for (std::size_t e=0; e<nelem; ++e)
365 : : {
366 : : // Reconstruct second-order dofs of volume-fractions in Taylor space
367 : : // using nodal-stencils, for a good interface-normal estimate
368 [ + - ]: 220760 : tk::recoLeastSqExtStencil( rdof, e, esup, inpoel, geoElem, P, vars );
369 : : }
370 : :
371 : : // 2.
372 [ + - ]: 1618 : tk::transform_P0P1(rdof, nelem, inpoel, coord, P, vars);
373 : 1618 : }
374 : :
375 : : //! Limit second-order solution, and primitive quantities separately
376 : : //! \param[in] geoFace Face geometry array
377 : : //! \param[in] fd Face connectivity and boundary conditions object
378 : : //! \param[in] esup Elements-surrounding-nodes connectivity
379 : : //! \param[in] inpoel Element-node connectivity
380 : : //! \param[in] coord Array of nodal coordinates
381 : : //! \param[in] srcFlag Whether the energy source was added
382 : : //! \param[in,out] U Solution vector at recent time step
383 : : //! \param[in,out] P Vector of primitives at recent time step
384 : 1618 : void limit( const tk::Fields& geoFace,
385 : : const inciter::FaceData& fd,
386 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
387 : : const std::vector< std::size_t >& inpoel,
388 : : const tk::UnsMesh::Coords& coord,
389 : : const std::vector< int >& srcFlag,
390 : : tk::Fields& U,
391 : : tk::Fields& P ) const
392 : : {
393 [ - + ][ - - ]: 1618 : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
394 : : "vector and primitive vector at recent time step incorrect" );
395 : :
396 : 1618 : const auto limiter = g_inputdeck.get< tag::limiter >();
397 : 1618 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
398 : : const auto& solidx = g_inputdeck.get<
399 : 1618 : tag::matidxmap, tag::solidx >();
400 : :
401 : : // limit vectors of conserved and primitive quantities
402 [ + - ]: 1618 : if (limiter == ctr::LimiterType::VERTEXBASEDP1)
403 : : {
404 : 1618 : VertexBasedMultiMat_FV( esup, inpoel, fd.Esuel().size()/4,
405 : : coord, srcFlag, solidx, U, P, nmat );
406 : 1618 : PositivityPreservingMultiMat_FV( inpoel, fd.Esuel().size()/4, nmat,
407 : 1618 : m_mat_blk, coord, geoFace, U, P );
408 : : }
409 [ - - ]: 0 : else if (limiter != ctr::LimiterType::NOLIMITER)
410 : : {
411 [ - - ][ - - ]: 0 : Throw("Limiter type not configured for multimat.");
[ - - ]
412 : : }
413 : 1618 : }
414 : :
415 : : //! Apply CPL to the conservative variable solution for this PDE system
416 : : //! \param[in] prim Array of primitive variables
417 : : //! \param[in] geoElem Element geometry array
418 : : //! \param[in] inpoel Element-node connectivity
419 : : //! \param[in] coord Array of nodal coordinates
420 : : //! \param[in,out] unk Array of conservative variables
421 : : //! \param[in] nielem Number of internal elements
422 : : //! \details This function applies CPL to obtain consistent dofs for
423 : : //! conservative quantities based on the limited primitive quantities.
424 : : //! See Pandare et al. (2023). On the Design of Stable,
425 : : //! Consistent, and Conservative High-Order Methods for Multi-Material
426 : : //! Hydrodynamics. J Comp Phys, 112313.
427 : : void CPL( const tk::Fields& prim,
428 : : const tk::Fields& geoElem,
429 : : const std::vector< std::size_t >& inpoel,
430 : : const tk::UnsMesh::Coords& coord,
431 : : tk::Fields& unk,
432 : : std::size_t nielem ) const
433 : : {
434 : : [[maybe_unused]] const auto rdof = g_inputdeck.get< tag::rdof >();
435 : : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
436 : :
437 : : Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
438 : : "vector and primitive vector at recent time step incorrect" );
439 : : Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
440 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
441 : : Assert( prim.nprop() == rdof*nprim(), "Number of components in vector of "
442 : : "primitive quantities must equal "+ std::to_string(rdof*nprim()) );
443 : :
444 : : correctLimConservMultiMat(nielem, m_mat_blk, nmat, inpoel,
445 : : coord, geoElem, prim, unk);
446 : : }
447 : :
448 : : //! Compute right hand side
449 : : //! \param[in] t Physical time
450 : : //! \param[in] geoFace Face geometry array
451 : : //! \param[in] geoElem Element geometry array
452 : : //! \param[in] fd Face connectivity and boundary conditions object
453 : : //! \param[in] inpoel Element-node connectivity
454 : : //! \param[in] coord Array of nodal coordinates
455 : : //! \param[in] elemblkid Element ids associated with mesh block ids where
456 : : //! user ICs are set
457 : : //! \param[in] U Solution vector at recent time step
458 : : //! \param[in] P Primitive vector at recent time step
459 : : //! \param[in,out] R Right-hand side vector computed
460 : : //! \param[in,out] srcFlag Whether the energy source was added
461 : 1618 : void rhs( tk::real t,
462 : : const tk::Fields& geoFace,
463 : : const tk::Fields& geoElem,
464 : : const inciter::FaceData& fd,
465 : : const std::vector< std::size_t >& inpoel,
466 : : const tk::UnsMesh::Coords& coord,
467 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
468 : : elemblkid,
469 : : const tk::Fields& U,
470 : : const tk::Fields& P,
471 : : tk::Fields& R,
472 : : std::vector< int >& srcFlag ) const
473 : : {
474 : 1618 : const auto rdof = g_inputdeck.get< tag::rdof >();
475 : 1618 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
476 : 1618 : const auto intsharp =
477 : 1618 : g_inputdeck.get< tag::multimat, tag::intsharp >();
478 : :
479 : 1618 : const auto nelem = fd.Esuel().size()/4;
480 : :
481 [ - + ][ - - ]: 1618 : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
482 : : "vector and primitive vector at recent time step incorrect" );
483 [ - + ][ - - ]: 1618 : Assert( U.nunk() == R.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
484 : : "vector and right-hand side at recent time step incorrect" );
485 [ - + ][ - - ]: 1618 : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
486 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
487 [ - + ][ - - ]: 1618 : Assert( P.nprop() == rdof*nprim(), "Number of components in primitive "
[ - - ][ - - ]
[ - - ]
488 : : "vector must equal "+ std::to_string(rdof*nprim()) );
489 [ - + ][ - - ]: 1618 : Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
[ - - ][ - - ]
490 : : "Mismatch in inpofa size" );
491 : :
492 : : // set rhs to zero
493 [ + - ]: 1618 : R.fill(0.0);
494 : :
495 : : // configure a no-op lambda for prescribed velocity
496 : 551116 : auto velfn = []( ncomp_t, tk::real, tk::real, tk::real, tk::real ){
497 : 551116 : return tk::VelFn::result_type(); };
498 : :
499 : : // compute internal surface flux (including non-conservative) integrals
500 [ + - ]: 3236 : tk::surfIntFV( nmat, m_mat_blk, t, rdof, inpoel,
501 [ + - ]: 1618 : coord, fd, geoFace, geoElem, m_riemann, velfn, U, P,
502 : : srcFlag, R, intsharp );
503 : :
504 : : // compute boundary surface flux (including non-conservative) integrals
505 [ + + ]: 11326 : for (const auto& b : m_bc)
506 : 9708 : tk::bndSurfIntFV( nmat, m_mat_blk, rdof, b.first,
507 [ + - ]: 9708 : fd, geoFace, geoElem, inpoel, coord, t, m_riemann,
508 [ + - ]: 9708 : velfn, b.second, U, P, srcFlag, R, intsharp );
509 : :
510 : : // compute optional source term
511 [ + - ][ + - ]: 1618 : tk::srcIntFV( m_mat_blk, t, fd.Esuel().size()/4,
512 : : geoElem, Problem::src, R, nmat );
513 : :
514 : : // compute finite pressure relaxation terms
515 [ + + ]: 1618 : if (g_inputdeck.get< tag::multimat, tag::prelax >())
516 : : {
517 : 1568 : const auto ct = g_inputdeck.get< tag::multimat,
518 : 1568 : tag::prelax_timescale >();
519 [ + - ]: 1568 : tk::pressureRelaxationIntFV( nmat, m_mat_blk, rdof,
520 : : nelem, inpoel, coord, geoElem, U, P, ct,
521 : : R );
522 : : }
523 : :
524 : : // compute external (energy) sources
525 [ + - ]: 1618 : m_physics.physSrc(nmat, t, geoElem, elemblkid, R, srcFlag);
526 : 1618 : }
527 : :
528 : : //! Compute the minimum time step size
529 : : // //! \param[in] fd Face connectivity and boundary conditions object
530 : : // //! \param[in] geoFace Face geometry array
531 : : //! \param[in] geoElem Element geometry array
532 : : //! \param[in] U Solution vector at recent time step
533 : : //! \param[in] P Vector of primitive quantities at recent time step
534 : : //! \param[in] nielem Number of internal elements
535 : : //! \param[in] srcFlag Whether the energy source was added
536 : : //! \param[in,out] local_dte Time step size for each element (for local
537 : : //! time stepping)
538 : : //! \return Minimum time step size
539 : : //! \details The allowable dt is calculated by looking at the maximum
540 : : //! wave-speed in elements surrounding each face, times the area of that
541 : : //! face. Once the maximum of this quantity over the mesh is determined,
542 : : //! the volume of each cell is divided by this quantity. A minimum of this
543 : : //! ratio is found over the entire mesh, which gives the allowable dt.
544 : 225 : tk::real dt( const inciter::FaceData& /*fd*/,
545 : : const tk::Fields& /*geoFace*/,
546 : : const tk::Fields& geoElem,
547 : : const tk::Fields& U,
548 : : const tk::Fields& P,
549 : : const std::size_t nielem,
550 : : const std::vector< int >& srcFlag,
551 : : std::vector< tk::real >& local_dte ) const
552 : : {
553 : 225 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
554 : :
555 : : // obtain dt restrictions from all physics
556 [ + - ]: 225 : auto dt_e = timeStepSizeMultiMatFV(m_mat_blk, geoElem, nielem, nmat, U,
557 : : P, local_dte);
558 [ + - ]: 225 : auto dt_p = m_physics.dtRestriction(geoElem, nielem, srcFlag);
559 : :
560 : 225 : return std::min(dt_e, dt_p);
561 : : }
562 : :
563 : : //! Extract the velocity field at cell nodes. Currently unused.
564 : : //! \param[in] U Solution vector at recent time step
565 : : //! \param[in] N Element node indices
566 : : //! \return Array of the four values of the velocity field
567 : : std::array< std::array< tk::real, 4 >, 3 >
568 : : velocity( const tk::Fields& U,
569 : : const std::array< std::vector< tk::real >, 3 >&,
570 : : const std::array< std::size_t, 4 >& N ) const
571 : : {
572 : : const auto rdof = g_inputdeck.get< tag::rdof >();
573 : : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
574 : :
575 : : std::array< std::array< tk::real, 4 >, 3 > v;
576 : : v[0] = U.extract( momentumDofIdx(nmat, 0, rdof, 0), N );
577 : : v[1] = U.extract( momentumDofIdx(nmat, 1, rdof, 0), N );
578 : : v[2] = U.extract( momentumDofIdx(nmat, 2, rdof, 0), N );
579 : :
580 : : std::vector< std::array< tk::real, 4 > > ar;
581 : : ar.resize(nmat);
582 : : for (std::size_t k=0; k<nmat; ++k)
583 : : ar[k] = U.extract( densityDofIdx(nmat, k, rdof, 0), N );
584 : :
585 : : std::array< tk::real, 4 > r{{ 0.0, 0.0, 0.0, 0.0 }};
586 : : for (std::size_t i=0; i<r.size(); ++i) {
587 : : for (std::size_t k=0; k<nmat; ++k)
588 : : r[i] += ar[k][i];
589 : : }
590 : :
591 : : std::transform( r.begin(), r.end(), v[0].begin(), v[0].begin(),
592 : : []( tk::real s, tk::real& d ){ return d /= s; } );
593 : : std::transform( r.begin(), r.end(), v[1].begin(), v[1].begin(),
594 : : []( tk::real s, tk::real& d ){ return d /= s; } );
595 : : std::transform( r.begin(), r.end(), v[2].begin(), v[2].begin(),
596 : : []( tk::real s, tk::real& d ){ return d /= s; } );
597 : : return v;
598 : : }
599 : :
600 : : //! Return a map that associates user-specified strings to functions
601 : : //! \return Map that associates user-specified strings to functions that
602 : : //! compute relevant quantities to be output to file
603 : 52 : std::map< std::string, tk::GetVarFn > OutVarFn() const
604 : 52 : { return MultiMatOutVarFn(); }
605 : :
606 : : //! Return analytic field names to be output to file
607 : : //! \return Vector of strings labelling analytic fields output in file
608 : 0 : std::vector< std::string > analyticFieldNames() const {
609 : 0 : auto nmat = g_inputdeck.get< eq, tag::nmat >();
610 : :
611 : 0 : return MultiMatFieldNames(nmat);
612 : : }
613 : :
614 : : //! Return surface field names to be output to file
615 : : //! \return Vector of strings labelling surface fields output in file
616 : 26 : std::vector< std::string > surfNames() const
617 : 26 : { return MultiMatSurfNames(); }
618 : :
619 : : //! Return time history field names to be output to file
620 : : //! \return Vector of strings labelling time history fields output in file
621 : 0 : std::vector< std::string > histNames() const {
622 : 0 : return MultiMatHistNames();
623 : : }
624 : :
625 : : //! Return surface field output going to file
626 : : std::vector< std::vector< tk::real > >
627 : 26 : surfOutput( const inciter::FaceData& fd,
628 : : const tk::Fields& U,
629 : : const tk::Fields& P ) const
630 : : {
631 : 26 : const auto rdof = g_inputdeck.get< tag::rdof >();
632 : 26 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
633 : :
634 : 26 : return MultiMatSurfOutput( nmat, rdof, fd, U, P );
635 : : }
636 : :
637 : : //! Return time history field output evaluated at time history points
638 : : //! \param[in] h History point data
639 : : //! \param[in] inpoel Element-node connectivity
640 : : //! \param[in] coord Array of nodal coordinates
641 : : //! \param[in] U Array of unknowns
642 : : //! \param[in] P Array of primitive quantities
643 : : //! \return Vector of time history output of bulk flow quantities (density,
644 : : //! velocity, total energy, pressure, and volume fraction) evaluated at
645 : : //! time history points
646 : : std::vector< std::vector< tk::real > >
647 : 0 : histOutput( const std::vector< HistData >& h,
648 : : const std::vector< std::size_t >& inpoel,
649 : : const tk::UnsMesh::Coords& coord,
650 : : const tk::Fields& U,
651 : : const tk::Fields& P ) const
652 : : {
653 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
654 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
655 : :
656 : 0 : const auto& x = coord[0];
657 : 0 : const auto& y = coord[1];
658 : 0 : const auto& z = coord[2];
659 : :
660 [ - - ]: 0 : std::vector< std::vector< tk::real > > Up(h.size());
661 : :
662 : 0 : std::size_t j = 0;
663 [ - - ]: 0 : for (const auto& p : h) {
664 : 0 : auto e = p.get< tag::elem >();
665 : 0 : auto chp = p.get< tag::coord >();
666 : :
667 : : // Evaluate inverse Jacobian
668 : 0 : std::array< std::array< tk::real, 3>, 4 > cp{{
669 : : {{ x[inpoel[4*e ]], y[inpoel[4*e ]], z[inpoel[4*e ]] }},
670 : 0 : {{ x[inpoel[4*e+1]], y[inpoel[4*e+1]], z[inpoel[4*e+1]] }},
671 : 0 : {{ x[inpoel[4*e+2]], y[inpoel[4*e+2]], z[inpoel[4*e+2]] }},
672 : 0 : {{ x[inpoel[4*e+3]], y[inpoel[4*e+3]], z[inpoel[4*e+3]] }} }};
673 : 0 : auto J = tk::inverseJacobian( cp[0], cp[1], cp[2], cp[3] );
674 : :
675 : : // evaluate solution at history-point
676 : 0 : std::array< tk::real, 3 > dc{{chp[0]-cp[0][0], chp[1]-cp[0][1],
677 : 0 : chp[2]-cp[0][2]}};
678 [ - - ]: 0 : auto B = tk::eval_basis(rdof, tk::dot(J[0],dc), tk::dot(J[1],dc),
679 : 0 : tk::dot(J[2],dc));
680 [ - - ]: 0 : auto uhp = eval_state(m_ncomp, rdof, rdof, e, U, B);
681 [ - - ]: 0 : auto php = eval_state(nprim(), rdof, rdof, e, P, B);
682 : :
683 : : // store solution in history output vector
684 [ - - ]: 0 : Up[j].resize(6+nmat, 0.0);
685 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
686 : 0 : Up[j][0] += uhp[densityIdx(nmat,k)];
687 : 0 : Up[j][4] += uhp[energyIdx(nmat,k)];
688 : 0 : Up[j][5] += php[pressureIdx(nmat,k)];
689 : 0 : Up[j][6+k] = uhp[volfracIdx(nmat,k)];
690 : : }
691 : 0 : Up[j][1] = php[velocityIdx(nmat,0)];
692 : 0 : Up[j][2] = php[velocityIdx(nmat,1)];
693 : 0 : Up[j][3] = php[velocityIdx(nmat,2)];
694 : 0 : ++j;
695 : : }
696 : :
697 : 0 : return Up;
698 : : }
699 : :
700 : : //! Return names of integral variables to be output to diagnostics file
701 : : //! \return Vector of strings labelling integral variables output
702 : 20 : std::vector< std::string > names() const
703 : : {
704 : 20 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
705 : 20 : return MultiMatDiagNames(nmat);
706 : : }
707 : :
708 : : //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
709 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
710 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
711 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
712 : : //! \param[in] t Physical time at which to evaluate the analytic solution
713 : : //! \return Vector of analytic solution at given location and time
714 : : std::vector< tk::real >
715 : 0 : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
716 : 0 : { return Problem::analyticSolution( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
717 : :
718 : : //! Return analytic solution for conserved variables
719 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
720 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
721 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
722 : : //! \param[in] t Physical time at which to evaluate the analytic solution
723 : : //! \return Vector of analytic solution at given location and time
724 : : std::vector< tk::real >
725 : 0 : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
726 : 0 : { return Problem::initialize( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
727 : :
728 : : //! Return cell-averaged specific total energy for an element
729 : : //! \param[in] e Element id for which total energy is required
730 : : //! \param[in] unk Vector of conserved quantities
731 : : //! \return Cell-averaged specific total energy for given element
732 : 0 : tk::real sp_totalenergy(std::size_t e, const tk::Fields& unk) const
733 : : {
734 : 0 : const auto rdof = g_inputdeck.get< tag::rdof >();
735 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
736 : :
737 : 0 : tk::real sp_te(0.0);
738 : : // sum each material total energy
739 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
740 : 0 : sp_te += unk(e, energyDofIdx(nmat,k,rdof,0));
741 : : }
742 : 0 : return sp_te;
743 : : }
744 : :
745 : : //! Compute relevant sound speed for output
746 : : //! \param[in] nielem Number of internal elements
747 : : //! \param[in] U Solution vector at recent time step
748 : : //! \param[in] P Primitive vector at recent time step
749 : : //! \param[in,out] ss Sound speed vector
750 : 26 : void soundspeed(
751 : : std::size_t nielem,
752 : : const tk::Fields& U,
753 : : const tk::Fields& P,
754 : : std::vector< tk::real >& ss) const
755 : : {
756 [ - + ][ - - ]: 26 : Assert( ss.size() == nielem, "Size of sound speed vector incorrect " );
[ - - ][ - - ]
757 : :
758 : 26 : const auto ndof = g_inputdeck.get< tag::ndof >();
759 : 26 : const auto rdof = g_inputdeck.get< tag::rdof >();
760 : 26 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
761 : 26 : std::size_t ncomp = U.nprop()/rdof;
762 : 26 : std::size_t nprim = P.nprop()/rdof;
763 : :
764 [ + - ][ + - ]: 52 : std::vector< tk::real > ugp(ncomp, 0.0), pgp(nprim, 0.0);
765 : :
766 [ + + ]: 8602 : for (std::size_t e=0; e<nielem; ++e) {
767 : : // basis function at centroid
768 [ + - ]: 17152 : std::vector< tk::real > B(rdof, 0.0);
769 : 8576 : B[0] = 1.0;
770 : :
771 : : // get conserved quantities
772 [ + - ]: 8576 : ugp = eval_state(ncomp, rdof, ndof, e, U, B);
773 : : // get primitive quantities
774 [ + - ]: 8576 : pgp = eval_state(nprim, rdof, ndof, e, P, B);
775 : :
776 : : // acoustic speed (this should be consistent with time-step calculation)
777 : 8576 : ss[e] = 0.0;
778 [ + + ]: 25728 : for (std::size_t k=0; k<nmat; ++k)
779 : : {
780 [ + + ]: 17152 : if (ugp[volfracIdx(nmat, k)] > 1.0e-04) {
781 [ + - ]: 21326 : ss[e] = std::max( ss[e], m_mat_blk[k].compute< EOS::soundspeed >(
782 : : ugp[densityIdx(nmat, k)], pgp[pressureIdx(nmat, k)],
783 : 10663 : ugp[volfracIdx(nmat, k)], k ) );
784 : : }
785 : : }
786 : : }
787 : 26 : }
788 : :
789 : : private:
790 : : //! Physics policy
791 : : const Physics m_physics;
792 : : //! Number of components in this PDE system
793 : : const ncomp_t m_ncomp;
794 : : //! Riemann solver
795 : : tk::RiemannFluxFn m_riemann;
796 : : //! BC configuration
797 : : BCStateFn m_bc;
798 : : //! EOS material block
799 : : std::vector< EOS > m_mat_blk;
800 : :
801 : : //! Evaluate conservative part of physical flux function for this PDE system
802 : : //! \param[in] ncomp Number of scalar components in this PDE system
803 : : //! \param[in] ugp Numerical solution at the Gauss point at which to
804 : : //! evaluate the flux
805 : : //! \return Flux vectors for all components in this PDE system
806 : : //! \note The function signature must follow tk::FluxFn
807 : : static tk::FluxFn::result_type
808 : : flux( ncomp_t ncomp,
809 : : const std::vector< EOS >& mat_blk,
810 : : const std::vector< tk::real >& ugp,
811 : : const std::vector< std::array< tk::real, 3 > >& )
812 : : {
813 : : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
814 : :
815 : : return tk::fluxTerms(ncomp, nmat, mat_blk, ugp);
816 : : }
817 : :
818 : : //! \brief Boundary state function providing the left and right state of a
819 : : //! face at Dirichlet boundaries
820 : : //! \param[in] ncomp Number of scalar components in this PDE system
821 : : //! \param[in] ul Left (domain-internal) state
822 : : //! \param[in] x X-coordinate at which to compute the states
823 : : //! \param[in] y Y-coordinate at which to compute the states
824 : : //! \param[in] z Z-coordinate at which to compute the states
825 : : //! \param[in] t Physical time
826 : : //! \return Left and right states for all scalar components in this PDE
827 : : //! system
828 : : //! \note The function signature must follow tk::StateFn. For multimat, the
829 : : //! left or right state is the vector of conserved quantities, followed by
830 : : //! the vector of primitive quantities appended to it.
831 : : static tk::StateFn::result_type
832 : 0 : dirichlet( ncomp_t ncomp,
833 : : const std::vector< EOS >& mat_blk,
834 : : const std::vector< tk::real >& ul, tk::real x, tk::real y,
835 : : tk::real z, tk::real t, const std::array< tk::real, 3 >& )
836 : : {
837 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
838 : :
839 [ - - ]: 0 : auto ur = Problem::initialize( ncomp, mat_blk, x, y, z, t );
840 [ - - ][ - - ]: 0 : Assert( ur.size() == ncomp, "Incorrect size for boundary state vector" );
[ - - ][ - - ]
841 : :
842 [ - - ]: 0 : ur.resize(ul.size());
843 : :
844 : 0 : tk::real rho(0.0);
845 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
846 : 0 : rho += ur[densityIdx(nmat, k)];
847 : :
848 : : // get primitives in boundary state
849 : :
850 : : // velocity
851 : 0 : ur[ncomp+velocityIdx(nmat, 0)] = ur[momentumIdx(nmat, 0)] / rho;
852 : 0 : ur[ncomp+velocityIdx(nmat, 1)] = ur[momentumIdx(nmat, 1)] / rho;
853 : 0 : ur[ncomp+velocityIdx(nmat, 2)] = ur[momentumIdx(nmat, 2)] / rho;
854 : :
855 : : // material pressures
856 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
857 : : {
858 [ - - ]: 0 : auto gk = getDeformGrad(nmat, k, ur);
859 [ - - ]: 0 : ur[ncomp+pressureIdx(nmat, k)] = mat_blk[k].compute< EOS::pressure >(
860 : 0 : ur[densityIdx(nmat, k)], ur[ncomp+velocityIdx(nmat, 0)],
861 : 0 : ur[ncomp+velocityIdx(nmat, 1)], ur[ncomp+velocityIdx(nmat, 2)],
862 : 0 : ur[energyIdx(nmat, k)], ur[volfracIdx(nmat, k)], k, gk );
863 : : }
864 : :
865 [ - - ][ - - ]: 0 : Assert( ur.size() == ncomp+nmat+3, "Incorrect size for appended "
[ - - ][ - - ]
866 : : "boundary state vector" );
867 : :
868 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
869 : : }
870 : :
871 : : // Other boundary condition types that do not depend on "Problem" should be
872 : : // added in BCFunctions.hpp
873 : : };
874 : :
875 : : } // fv::
876 : :
877 : : } // inciter::
878 : :
879 : : #endif // FVMultiMat_h
|
