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