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