Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/MultiMat/BCFunctions.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 Functions specifying boundary conditions.
9 : : \details Functions that return boundary state when the interior state at
10 : : at the boundary location is provided.
11 : : */
12 : : // *****************************************************************************
13 : : #ifndef BCFunctions_h
14 : : #define BCFunctions_h
15 : :
16 : : #include "FunctionPrototypes.hpp"
17 : : #include "MiscMultiMatFns.hpp"
18 : :
19 : : namespace inciter {
20 : :
21 : : //! \brief Boundary state function providing the left and right state of a
22 : : //! face at symmetry boundaries
23 : : //! \param[in] ncomp Number of scalar components in this PDE system
24 : : //! \param[in] ul Left (domain-internal) state
25 : : //! \param[in] fn Unit face normal
26 : : //! \return Left and right states for all scalar components in this PDE
27 : : //! system
28 : : //! \note The function signature must follow tk::StateFn. For multimat, the
29 : : //! left or right state is the vector of conserved quantities, followed by
30 : : //! the vector of primitive quantities appended to it.
31 : : static tk::StateFn::result_type
32 : 1791280 : symmetry( ncomp_t ncomp,
33 : : const std::vector< EOS >&,
34 : : const std::vector< tk::real >& ul,
35 : : tk::real, tk::real, tk::real, tk::real,
36 : : const std::array< tk::real, 3 >& fn )
37 : : {
38 : 1791280 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
39 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
40 : :
41 : 1791280 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
42 : :
43 : : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
44 : : "internal state vector" );
45 : :
46 : : tk::real rho(0.0);
47 [ + + ]: 5997240 : for (std::size_t k=0; k<nmat; ++k)
48 : 4205960 : rho += ul[densityIdx(nmat, k)];
49 : :
50 : 1791280 : auto ur = ul;
51 : :
52 : : // Internal cell velocity components
53 : 1791280 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
54 : 1791280 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
55 : 1791280 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
56 : : // Normal component of velocity
57 : 1791280 : auto vnl = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
58 : : // Ghost state velocity components
59 : 1791280 : auto v1r = v1l - 2.0*vnl*fn[0];
60 : 1791280 : auto v2r = v2l - 2.0*vnl*fn[1];
61 : 1791280 : auto v3r = v3l - 2.0*vnl*fn[2];
62 : : // Boundary condition
63 [ + + ]: 5997240 : for (std::size_t k=0; k<nmat; ++k)
64 : : {
65 [ - + ]: 4205960 : ur[volfracIdx(nmat, k)] = ul[volfracIdx(nmat, k)];
66 [ - + ]: 4205960 : ur[densityIdx(nmat, k)] = ul[densityIdx(nmat, k)];
67 : 4205960 : ur[energyIdx(nmat, k)] = ul[energyIdx(nmat, k)];
68 [ - + ]: 4205960 : if (solidx[k] > 0) {
69 : : // Internal inverse deformation tensor
70 : : std::array< std::array< tk::real, 3 >, 3 > g;
71 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
72 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
73 : 0 : g[i][j] = ul[deformIdx(nmat,solidx[k],i,j)];
74 : : // Internal Cauchy stress tensor
75 : : std::array< std::array< tk::real, 3 >, 3 > s;
76 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
77 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
78 : 0 : s[i][j] = ul[ncomp+stressIdx(nmat,solidx[k],stressCmp[i][j])];
79 : : // Make reflection matrix
80 : : std::array< std::array< tk::real, 3 >, 3 >
81 : 0 : reflectionMat{{{1,0,0}, {0,1,0}, {0,0,1}}};
82 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
83 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
84 : 0 : reflectionMat[i][j] -= 2*fn[i]*fn[j];
85 : : // Reflect g
86 [ - - ]: 0 : g = tk::reflectTensor(g, reflectionMat);
87 : : // Reflect s
88 [ - - ]: 0 : s = tk::reflectTensor(s, reflectionMat);
89 : : // Copy g and s into ur
90 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
91 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
92 : 0 : ur[deformIdx(nmat,solidx[k],i,j)] = g[i][j];
93 : 0 : ur[ncomp+stressIdx(nmat,solidx[k],stressCmp[i][j])] = s[i][j];
94 : : }
95 : : }
96 : : }
97 : 1791280 : ur[momentumIdx(nmat, 0)] = rho * v1r;
98 : 1791280 : ur[momentumIdx(nmat, 1)] = rho * v2r;
99 : 1791280 : ur[momentumIdx(nmat, 2)] = rho * v3r;
100 : :
101 : : // Internal cell primitive quantities using the separately reconstructed
102 : : // primitive quantities. This is used to get ghost state for primitive
103 : : // quantities
104 : :
105 : : // velocity
106 : 1791280 : ur[ncomp+velocityIdx(nmat, 0)] = v1r;
107 : 1791280 : ur[ncomp+velocityIdx(nmat, 1)] = v2r;
108 : 1791280 : ur[ncomp+velocityIdx(nmat, 2)] = v3r;
109 : : // material pressures
110 [ + + ]: 5997240 : for (std::size_t k=0; k<nmat; ++k)
111 : 4205960 : ur[ncomp+pressureIdx(nmat, k)] = ul[ncomp+pressureIdx(nmat, k)];
112 : :
113 : : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
114 : : "boundary state vector" );
115 : :
116 [ + - ]: 1791280 : return {{ std::move(ul), std::move(ur) }};
117 : : }
118 : :
119 : : //! \brief Boundary state function providing the left and right state of a
120 : : //! face at inlet boundaries
121 : : //! \param[in] ncomp Number of scalar components in this PDE system
122 : : //! \param[in] ul Left (domain-internal) state
123 : : //! \param[in] fn Unit face normal
124 : : //! \return Left and right states for all scalar components in this PDE
125 : : //! system
126 : : //! \details The inlet boundary condition specifies a velocity at a
127 : : //! sideset and assumes a zero bulk pressure and density gradient
128 : : //! \note The function signature must follow tk::StateFn
129 : : static tk::StateFn::result_type
130 : 0 : inlet( ncomp_t ncomp,
131 : : const std::vector< EOS >& mat_blk,
132 : : const std::vector< tk::real >& ul,
133 : : tk::real, tk::real, tk::real, tk::real,
134 : : const std::array< tk::real, 3 >& fn )
135 : : {
136 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
137 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
138 : 0 : auto& inbc = g_inputdeck.get< tag::bc >()[0].get< tag::inlet >();
139 : :
140 : : // inlet velocity and material
141 : 0 : auto u_in = inbc[0].get< tag::velocity >();
142 : 0 : auto mat_in = inbc[0].get< tag::materialid >() - 1;
143 : 0 : auto p_in = inbc[0].get< tag::pressure >();
144 : 0 : auto t_in = inbc[0].get< tag::temperature >();
145 : :
146 [ - - ]: 0 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
147 : :
148 : : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
149 : : "internal state vector" );
150 : :
151 [ - - ]: 0 : auto ur = ul;
152 : :
153 : : // External cell velocity, such that velocity = v_in at face
154 : 0 : auto v1r = u_in[0];
155 : 0 : auto v2r = u_in[1];
156 : 0 : auto v3r = u_in[2];
157 : :
158 : : // Normal inlet velocity
159 : 0 : auto vn = u_in[0]*fn[0] + u_in[1]*fn[1] + u_in[2]*fn[2];
160 : :
161 : : // Acoustic speed
162 : 0 : tk::real a(0.0);
163 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
164 [ - - ]: 0 : if (ul[volfracIdx(nmat, k)] > 1.0e-04)
165 [ - - ][ - - ]: 0 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
166 [ - - ]: 0 : ul[densityIdx(nmat, k)], ul[ncomp+pressureIdx(nmat, k)],
167 : : ul[volfracIdx(nmat, k)], k ) );
168 : :
169 : : // Mach number
170 : 0 : auto Ma = vn / a;
171 : :
172 : : tk::real alphamin = 1e-12;
173 : 0 : tk::real pk(0.0);
174 : : tk::real rho(0.0);
175 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
176 [ - - ]: 0 : if (k == mat_in)
177 : 0 : ur[volfracIdx(nmat,k)] = 1.0 -
178 : 0 : (static_cast< tk::real >(nmat-1))*alphamin;
179 : : else
180 : 0 : ur[volfracIdx(nmat,k)] = alphamin;
181 : :
182 : : // Material pressure, which, for supersonic inflow, is the exterior
183 : : // pressure and the interior pressure for subsonic
184 [ - - ]: 0 : if(Ma <= -1)
185 : 0 : pk = p_in;
186 : : else
187 : 0 : pk = ul[ncomp+pressureIdx(nmat,k)]/ul[volfracIdx(nmat,k)];
188 [ - - ]: 0 : auto rhok = mat_blk[k].compute< EOS::density >(pk, t_in);
189 : :
190 [ - - ]: 0 : ur[ncomp+pressureIdx(nmat, k)] = ur[volfracIdx(nmat,k)] * pk;
191 [ - - ]: 0 : ur[densityIdx(nmat,k)] = ur[volfracIdx(nmat,k)] * rhok;
192 : 0 : ur[energyIdx(nmat,k)] =
193 [ - - ]: 0 : mat_blk[k].compute< EOS::totalenergy >(ur[volfracIdx(nmat,k)]*rhok,
194 [ - - ][ - - ]: 0 : v1r, v2r, v3r, ur[volfracIdx(nmat,k)]*pk, ur[volfracIdx(nmat,k)]);
195 : :
196 : : // bulk density
197 : 0 : rho += ur[densityIdx(nmat,k)];
198 : : }
199 : :
200 [ - - ]: 0 : ur[momentumIdx(nmat, 0)] = rho * v1r;
201 : 0 : ur[momentumIdx(nmat, 1)] = rho * v2r;
202 : 0 : ur[momentumIdx(nmat, 2)] = rho * v3r;
203 : :
204 : : // velocity
205 : 0 : ur[ncomp+velocityIdx(nmat, 0)] = v1r;
206 : 0 : ur[ncomp+velocityIdx(nmat, 1)] = v2r;
207 [ - - ]: 0 : ur[ncomp+velocityIdx(nmat, 2)] = v3r;
208 : :
209 : : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
210 : : "boundary state vector" );
211 : :
212 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
213 : : }
214 : :
215 : : //! \brief Boundary state function providing the left and right state of a
216 : : //! face at farfield boundaries
217 : : //! \param[in] ncomp Number of scalar components in this PDE system
218 : : //! \param[in] ul Left (domain-internal) state
219 : : //! \param[in] fn Unit face normal
220 : : //! \return Left and right states for all scalar components in this PDE
221 : : //! system
222 : : //! \details The farfield boudary calculation, implemented here, is
223 : : //! based on the characteristic theory of hyperbolic systems.
224 : : //! \note The function signature must follow tk::StateFn
225 : : static tk::StateFn::result_type
226 : 0 : farfield( ncomp_t ncomp,
227 : : const std::vector< EOS >& mat_blk,
228 : : const std::vector< tk::real >& ul,
229 : : tk::real, tk::real, tk::real, tk::real,
230 : : const std::array< tk::real, 3 >& fn )
231 : : {
232 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
233 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
234 : :
235 : : // Farfield primitive quantities
236 : : auto fp =
237 : 0 : g_inputdeck.get< tag::bc >()[0].get< tag::pressure >();
238 : : auto ft =
239 : 0 : g_inputdeck.get< tag::bc >()[0].get< tag::temperature >();
240 : : auto fu =
241 : 0 : g_inputdeck.get< tag::bc >()[0].get< tag::velocity >();
242 : : auto fmat =
243 : 0 : g_inputdeck.get< tag::bc >()[0].get< tag::materialid >() - 1;
244 : :
245 [ - - ]: 0 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
246 : :
247 : : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
248 : : "internal state vector" );
249 : :
250 [ - - ]: 0 : auto ur = ul;
251 : :
252 : : // Internal cell velocity components
253 : 0 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
254 : 0 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
255 : 0 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
256 : :
257 : : // Normal velocity
258 : 0 : auto vn = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
259 : :
260 : : // Acoustic speed
261 : 0 : tk::real a(0.0);
262 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
263 [ - - ]: 0 : if (ul[volfracIdx(nmat, k)] > 1.0e-04)
264 [ - - ][ - - ]: 0 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
265 [ - - ]: 0 : ul[densityIdx(nmat, k)], ul[ncomp+pressureIdx(nmat, k)],
266 : : ul[volfracIdx(nmat, k)], k ) );
267 : :
268 : : // Mach number
269 : 0 : auto Ma = vn / a;
270 : :
271 : : tk::real alphamin = 1e-12;
272 : :
273 [ - - ]: 0 : if (Ma <= -1) { // Supersonic inflow
274 : : // For supersonic inflow, all the characteristics are from outside.
275 : : // Therefore, we calculate the ghost cell state using the primitive
276 : : // variables from outside.
277 : : tk::real rho(0.0);
278 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
279 [ - - ]: 0 : if (k == fmat)
280 : 0 : ur[volfracIdx(nmat,k)] = 1.0 -
281 : 0 : (static_cast< tk::real >(nmat-1))*alphamin;
282 : : else
283 : 0 : ur[volfracIdx(nmat,k)] = alphamin;
284 [ - - ]: 0 : auto rhok = mat_blk[k].compute< EOS::density >(fp, ft);
285 [ - - ]: 0 : ur[densityIdx(nmat,k)] = ur[volfracIdx(nmat,k)] * rhok;
286 : 0 : ur[energyIdx(nmat,k)] =
287 [ - - ]: 0 : mat_blk[k].compute< EOS::totalenergy >(ur[volfracIdx(nmat,k)]*rhok,
288 [ - - ]: 0 : fu[0], fu[1], fu[2], ur[volfracIdx(nmat,k)]*fp, ur[volfracIdx(nmat,k)]);
289 : :
290 : : // material pressures
291 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ur[volfracIdx(nmat, k)] * fp;
292 : :
293 : 0 : rho += ur[densityIdx(nmat,k)];
294 : : }
295 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
296 : 0 : ur[momentumIdx(nmat,i)] = rho * fu[i];
297 : 0 : ur[ncomp+velocityIdx(nmat, i)] = fu[i];
298 : : }
299 : :
300 [ - - ][ - - ]: 0 : } else if (Ma > -1 && Ma < 0) { // Subsonic inflow
301 : : // For subsonic inflow, there is 1 outgoing characteristic and 4
302 : : // incoming characteristics. Therefore, we calculate the ghost cell state
303 : : // by taking pressure from the internal cell and other quantities from
304 : : // the outside.
305 : : tk::real rho(0.0);
306 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
307 [ - - ]: 0 : if (k == fmat)
308 : 0 : ur[volfracIdx(nmat,k)] = 1.0 -
309 : 0 : (static_cast< tk::real >(nmat-1))*alphamin;
310 : : else
311 : 0 : ur[volfracIdx(nmat,k)] = alphamin;
312 [ - - ]: 0 : auto p = ul[ncomp+pressureIdx(nmat,k)] / ul[volfracIdx(nmat,k)];
313 [ - - ]: 0 : auto rhok = mat_blk[k].compute< EOS::density >(p, ft);
314 [ - - ]: 0 : ur[densityIdx(nmat,k)] = ur[volfracIdx(nmat,k)] * rhok;
315 : 0 : ur[energyIdx(nmat,k)] =
316 [ - - ]: 0 : mat_blk[k].compute< EOS::totalenergy >(ur[volfracIdx(nmat,k)]*rhok,
317 [ - - ]: 0 : fu[0], fu[1], fu[2], ur[volfracIdx(nmat,k)]*p, ur[volfracIdx(nmat,k)]);
318 : :
319 : : // material pressures
320 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ur[volfracIdx(nmat, k)] * p;
321 : :
322 : 0 : rho += ur[densityIdx(nmat,k)];
323 : : }
324 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
325 : 0 : ur[momentumIdx(nmat,i)] = rho * fu[i];
326 : 0 : ur[ncomp+velocityIdx(nmat, i)] = fu[i];
327 : : }
328 : :
329 [ - - ][ - - ]: 0 : } else if (Ma >= 0 && Ma < 1) { // Subsonic outflow
330 : : // For subsonic outflow, there is 1 incoming characteristic and 4
331 : : // outgoing characteristics. Therefore, we calculate the ghost cell state
332 : : // by taking pressure from the outside and other quantities from the
333 : : // internal cell.
334 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
335 : 0 : ur[energyIdx(nmat, k)] =
336 [ - - ]: 0 : mat_blk[k].compute< EOS::totalenergy >(
337 [ - - ]: 0 : ur[densityIdx(nmat, k)], v1l, v2l, v3l, ul[volfracIdx(nmat, k)]*fp,
338 [ - - ]: 0 : ul[volfracIdx(nmat, k)] );
339 : :
340 : : // material pressures
341 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ul[volfracIdx(nmat, k)] * fp;
342 : : }
343 : : }
344 : : // Otherwise, for supersonic outflow, all the characteristics are from
345 : : // internal cell. Therefore, we calculate the ghost cell state using the
346 : : // conservative variables from internal cell (which is what ur is
347 : : // initialized to).
348 : :
349 : : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
350 : : "boundary state vector" );
351 : :
352 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
353 : : }
354 : :
355 : : //! \brief Boundary state function providing the left and right state of a
356 : : //! face at extrapolation boundaries
357 : : //! \param[in] ul Left (domain-internal) state
358 : : //! \return Left and right states for all scalar components in this PDE
359 : : //! system
360 : : //! \note The function signature must follow tk::StateFn. For multimat, the
361 : : //! left or right state is the vector of conserved quantities, followed by
362 : : //! the vector of primitive quantities appended to it.
363 : : static tk::StateFn::result_type
364 : 107476 : extrapolate( ncomp_t,
365 : : const std::vector< EOS >&,
366 : : const std::vector< tk::real >& ul,
367 : : tk::real, tk::real, tk::real, tk::real,
368 : : const std::array< tk::real, 3 >& )
369 : : {
370 : 107476 : return {{ ul, ul }};
371 : : }
372 : :
373 : : //! \brief Boundary state function providing the left and right state of a
374 : : //! face at no-slip wall boundaries
375 : : //! \param[in] ncomp Number of scalar components in this PDE system
376 : : //! \param[in] ul Left (domain-internal) state
377 : : //! \param[in] fn Unit face normal
378 : : //! \return Left and right states for all scalar components in this PDE
379 : : //! system
380 : : //! \note The function signature must follow tk::StateFn. For multimat, the
381 : : //! left or right state is the vector of conserved quantities, followed by
382 : : //! the vector of primitive quantities appended to it.
383 : : static tk::StateFn::result_type
384 : 0 : noslipwall( ncomp_t ncomp,
385 : : const std::vector< EOS >&,
386 : : const std::vector< tk::real >& ul,
387 : : tk::real, tk::real, tk::real, tk::real,
388 : : const std::array< tk::real, 3 >& fn )
389 : : {
390 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
391 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
392 : :
393 : 0 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
394 : :
395 : : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
396 : : "internal state vector" );
397 : :
398 : : tk::real rho(0.0);
399 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
400 : 0 : rho += ul[densityIdx(nmat, k)];
401 : :
402 : 0 : auto ur = ul;
403 : :
404 : : // Internal cell velocity components
405 : 0 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
406 : 0 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
407 : 0 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
408 : : // Ghost state velocity components
409 : 0 : auto v1r = -v1l;
410 : 0 : auto v2r = -v2l;
411 : 0 : auto v3r = -v3l;
412 : : // Boundary condition
413 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
414 : : {
415 [ - - ]: 0 : ur[volfracIdx(nmat, k)] = ul[volfracIdx(nmat, k)];
416 [ - - ]: 0 : ur[densityIdx(nmat, k)] = ul[densityIdx(nmat, k)];
417 : 0 : ur[energyIdx(nmat, k)] = ul[energyIdx(nmat, k)];
418 [ - - ]: 0 : if (solidx[k] > 0) {
419 : : // Internal inverse deformation tensor
420 : : std::array< std::array< tk::real, 3 >, 3 > g;
421 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
422 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
423 : 0 : g[i][j] = ul[deformIdx(nmat,solidx[k],i,j)];
424 : : // Internal Cauchy stress tensor
425 : : std::array< std::array< tk::real, 3 >, 3 > s;
426 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
427 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
428 : 0 : s[i][j] = ul[ncomp+stressIdx(nmat,solidx[k],stressCmp[i][j])];
429 : : // Make reflection matrix
430 : : std::array< std::array< tk::real, 3 >, 3 >
431 : 0 : reflectionMat{{{1,0,0}, {0,1,0}, {0,0,1}}};
432 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
433 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
434 : 0 : reflectionMat[i][j] -= 2*fn[i]*fn[j];
435 : : // Reflect g
436 [ - - ]: 0 : g = tk::reflectTensor(g, reflectionMat);
437 : : // Reflect s
438 [ - - ]: 0 : s = tk::reflectTensor(s, reflectionMat);
439 : : // Copy g and s into ur
440 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
441 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
442 : 0 : ur[deformIdx(nmat,solidx[k],i,j)] = g[i][j];
443 : 0 : ur[ncomp+stressIdx(nmat,solidx[k],stressCmp[i][j])] = s[i][j];
444 : : }
445 : : }
446 : : }
447 : 0 : ur[momentumIdx(nmat, 0)] = rho * v1r;
448 : 0 : ur[momentumIdx(nmat, 1)] = rho * v2r;
449 : 0 : ur[momentumIdx(nmat, 2)] = rho * v3r;
450 : :
451 : : // Internal cell primitive quantities using the separately reconstructed
452 : : // primitive quantities. This is used to get ghost state for primitive
453 : : // quantities
454 : :
455 : : // velocity
456 : 0 : ur[ncomp+velocityIdx(nmat, 0)] = v1r;
457 : 0 : ur[ncomp+velocityIdx(nmat, 1)] = v2r;
458 : 0 : ur[ncomp+velocityIdx(nmat, 2)] = v3r;
459 : : // material pressures
460 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
461 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ul[ncomp+pressureIdx(nmat, k)];
462 : :
463 : : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
464 : : "boundary state vector" );
465 : :
466 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
467 : : }
468 : :
469 : : //! \brief Boundary state function providing the left and right state of a
470 : : //! face at back pressure boundaries
471 : : //! \param[in] ncomp Number of scalar components in this PDE system
472 : : //! \param[in] ul Left (domain-internal) state
473 : : //! \return Left and right states for all scalar components in this PDE
474 : : //! system
475 : : //! \note The function signature must follow tk::StateFn
476 : : static tk::StateFn::result_type
477 : 0 : back_pressure( ncomp_t ncomp,
478 : : const std::vector< EOS >& mat_blk,
479 : : const std::vector< tk::real >& ul,
480 : : tk::real, tk::real, tk::real, tk::real,
481 : : const std::array< tk::real, 3 >& )
482 : : {
483 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
484 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
485 : :
486 : : // Back pressure
487 : 0 : auto fbp = g_inputdeck.get< tag::bc >()[0].get< tag::back_pressure >().get<
488 : 0 : tag::pressure >();
489 : :
490 : 0 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
491 : :
492 : : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
493 : : "internal state vector" );
494 : :
495 : 0 : auto ur = ul;
496 : :
497 : : // Internal cell velocity components
498 : : std::array< tk::real, 3 > vl;
499 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
500 : 0 : vl[i] = ul[ncomp+velocityIdx(nmat, i)];
501 : :
502 : : // The ghost cell state is calculated using the back pressure and other
503 : : // quantities from the internal cell
504 : : auto rhor_b(0.0);
505 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
506 [ - - ]: 0 : auto T_k = mat_blk[k].compute< inciter::EOS::temperature >(
507 : : ul[densityIdx(nmat, k)], vl[0], vl[1], vl[2], ul[energyIdx(nmat, k)],
508 [ - - ]: 0 : ul[volfracIdx(nmat, k)] );
509 [ - - ]: 0 : auto rhor_k = mat_blk[k].compute< inciter::EOS::density >( fbp, T_k );
510 [ - - ]: 0 : ur[densityIdx(nmat, k)] = ul[volfracIdx(nmat, k)] * rhor_k;
511 : 0 : ur[energyIdx(nmat, k)] =
512 : 0 : mat_blk[k].compute< EOS::totalenergy >(
513 : 0 : ul[volfracIdx(nmat, k)]*rhor_k, vl[0], vl[1], vl[2],
514 [ - - ][ - - ]: 0 : ul[volfracIdx(nmat, k)]*fbp, ul[volfracIdx(nmat, k)] );
515 : 0 : rhor_b += ur[densityIdx(nmat, k)];
516 : :
517 : : // material pressures
518 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ul[volfracIdx(nmat, k)] * fbp;
519 : : }
520 : :
521 : : // bulk momentum
522 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
523 : 0 : ur[momentumIdx(nmat,i)] = rhor_b * vl[i];
524 : :
525 : : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
526 : : "boundary state vector" );
527 : :
528 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
529 : : }
530 : :
531 : : //----------------------------------------------------------------------------
532 : : // Boundary Gradient functions
533 : : //----------------------------------------------------------------------------
534 : :
535 : : //! \brief Boundary gradient function copying the left gradient to the right
536 : : //! gradient at a face
537 : : //! \param[in] dul Left (domain-internal) state
538 : : //! \return Left and right states for all scalar components in this PDE
539 : : //! system
540 : : //! \note The function signature must follow tk::StateFn. For multimat, the
541 : : //! left or right state is the vector of gradients of primitive quantities.
542 : : static tk::StateFn::result_type
543 : 0 : noOpGrad( ncomp_t,
544 : : const std::vector< EOS >&,
545 : : const std::vector< tk::real >& dul,
546 : : tk::real, tk::real, tk::real, tk::real,
547 : : const std::array< tk::real, 3 >& )
548 : : {
549 : 0 : return {{ dul, dul }};
550 : : }
551 : :
552 : : //! \brief Boundary gradient function for the symmetry boundary condition
553 : : //! \param[in] ncomp Number of variables whos gradients are needed
554 : : //! \param[in] dul Left (domain-internal) gradients
555 : : //! \return Left and right states for all scalar components in this PDE
556 : : //! system
557 : : //! \note The function signature must follow tk::StateFn. For multimat, the
558 : : //! left or right state is the vector of gradients of primitive quantities.
559 : : static tk::StateFn::result_type
560 : 0 : symmetryGrad( ncomp_t ncomp,
561 : : const std::vector< EOS >&,
562 : : const std::vector< tk::real >& dul,
563 : : tk::real, tk::real, tk::real, tk::real,
564 : : const std::array< tk::real, 3 >& )
565 : : {
566 : : Assert(dul.size() == 3*ncomp, "Incorrect size of boundary gradient vector");
567 : :
568 : 0 : auto dur = dul;
569 : :
570 [ - - ]: 0 : for (std::size_t i=0; i<3*ncomp; ++i)
571 : 0 : dur[i] = -dul[i];
572 : :
573 [ - - ]: 0 : return {{ std::move(dul), std::move(dur) }};
574 : : }
575 : :
576 : : //! \brief Boundary gradient function for zero gradient cells
577 : : //! \param[in] ncomp Number of variables whos gradients are needed
578 : : //! \param[in] dul Left (domain-internal) gradients
579 : : //! \return Left and right states for all scalar components in this PDE
580 : : //! system
581 : : //! \note The function signature must follow tk::StateFn. For multimat, the
582 : : //! left or right state is the vector of gradients of primitive quantities.
583 : : static tk::StateFn::result_type
584 : 0 : zeroGrad( ncomp_t ncomp,
585 : : const std::vector< EOS >&,
586 : : const std::vector< tk::real >& dul,
587 : : tk::real, tk::real, tk::real, tk::real,
588 : : const std::array< tk::real, 3 >& )
589 : : {
590 : : Assert(dul.size() == 3*ncomp, "Incorrect size of boundary gradient vector");
591 : :
592 : 0 : std::vector< tk::real > dur(3*ncomp, 0.0);
593 : :
594 [ - - ]: 0 : return {{ std::move(dul), std::move(dur) }};
595 : : }
596 : :
597 : : } // inciter::
598 : :
599 : : #endif // BCFunctions_h
|
