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 : 1983280 : 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 : 1983280 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
39 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
40 : :
41 : 1983280 : [[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 [ + + ]: 6573240 : for (std::size_t k=0; k<nmat; ++k)
48 : 4589960 : rho += ul[densityIdx(nmat, k)];
49 : :
50 : 1983280 : auto ur = ul;
51 : :
52 : : // Internal cell velocity components
53 : 1983280 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
54 : 1983280 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
55 : 1983280 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
56 : : // Normal component of velocity
57 : 1983280 : auto vnl = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
58 : : // Ghost state velocity components
59 : 1983280 : auto v1r = v1l - 2.0*vnl*fn[0];
60 : 1983280 : auto v2r = v2l - 2.0*vnl*fn[1];
61 : 1983280 : auto v3r = v3l - 2.0*vnl*fn[2];
62 : : // Boundary condition
63 [ + + ]: 6573240 : for (std::size_t k=0; k<nmat; ++k)
64 : : {
65 [ - + ]: 4589960 : ur[volfracIdx(nmat, k)] = ul[volfracIdx(nmat, k)];
66 [ - + ]: 4589960 : ur[densityIdx(nmat, k)] = ul[densityIdx(nmat, k)];
67 : 4589960 : ur[energyIdx(nmat, k)] = ul[energyIdx(nmat, k)];
68 [ - + ]: 4589960 : 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 : 1983280 : ur[momentumIdx(nmat, 0)] = rho * v1r;
98 : 1983280 : ur[momentumIdx(nmat, 1)] = rho * v2r;
99 : 1983280 : 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 : 1983280 : ur[ncomp+velocityIdx(nmat, 0)] = v1r;
107 : 1983280 : ur[ncomp+velocityIdx(nmat, 1)] = v2r;
108 : 1983280 : ur[ncomp+velocityIdx(nmat, 2)] = v3r;
109 : : // material pressures
110 [ + + ]: 6573240 : for (std::size_t k=0; k<nmat; ++k)
111 : 4589960 : 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 [ + - ]: 1983280 : 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 : 0 : tk::real alphamin = g_inputdeck.get< tag::multimat, tag::min_volumefrac >();
140 : :
141 : : // inlet velocity and material
142 : 0 : auto u_in = inbc[0].get< tag::velocity >();
143 : 0 : auto mat_in = inbc[0].get< tag::materialid >() - 1;
144 : 0 : auto p_in = inbc[0].get< tag::pressure >();
145 : 0 : auto t_in = inbc[0].get< tag::temperature >();
146 : :
147 [ - - ]: 0 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
148 : :
149 : : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
150 : : "internal state vector" );
151 : :
152 [ - - ]: 0 : auto ur = ul;
153 : :
154 : : // External cell velocity, such that velocity = v_in at face
155 : 0 : auto v1r = u_in[0];
156 : 0 : auto v2r = u_in[1];
157 : 0 : auto v3r = u_in[2];
158 : :
159 : : // Normal inlet velocity
160 : 0 : auto vn = u_in[0]*fn[0] + u_in[1]*fn[1] + u_in[2]*fn[2];
161 : :
162 : : // Acoustic speed
163 : 0 : tk::real a(0.0);
164 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
165 [ - - ]: 0 : if (ul[volfracIdx(nmat, k)] > 1.0e-04)
166 [ - - ][ - - ]: 0 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
167 [ - - ]: 0 : ul[densityIdx(nmat, k)], ul[ncomp+pressureIdx(nmat, k)],
168 : : ul[volfracIdx(nmat, k)], k ) );
169 : :
170 : : // Mach number
171 : 0 : auto Ma = vn / a;
172 : :
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 : 1500 : 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 : 1500 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
233 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
234 : 1500 : tk::real alphamin = g_inputdeck.get< tag::multimat, tag::min_volumefrac >();
235 : :
236 : : // Farfield primitive quantities
237 : : auto fp =
238 : 1500 : g_inputdeck.get< tag::bc >()[0].get< tag::pressure >();
239 : : auto ft =
240 : 1500 : g_inputdeck.get< tag::bc >()[0].get< tag::temperature >();
241 : : auto fu =
242 : 1500 : g_inputdeck.get< tag::bc >()[0].get< tag::velocity >();
243 : : auto fmat =
244 : 1500 : g_inputdeck.get< tag::bc >()[0].get< tag::materialid >() - 1;
245 : :
246 [ + - ]: 1500 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
247 : :
248 : : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
249 : : "internal state vector" );
250 : :
251 [ + - ]: 1500 : auto ur = ul;
252 : :
253 : : // Internal cell velocity components
254 : 1500 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
255 : 1500 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
256 : 1500 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
257 : :
258 : : // Normal velocity
259 : 1500 : auto vn = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
260 : :
261 : : // Acoustic speed
262 : 1500 : tk::real a(0.0);
263 [ + + ]: 4500 : for (std::size_t k=0; k<nmat; ++k)
264 [ + + ]: 3000 : if (ul[volfracIdx(nmat, k)] > 1.0e-04)
265 [ + - ][ + - ]: 3000 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
266 [ + - ]: 1500 : ul[densityIdx(nmat, k)], ul[ncomp+pressureIdx(nmat, k)],
267 : : ul[volfracIdx(nmat, k)], k ) );
268 : :
269 : : // Mach number
270 : 1500 : auto Ma = vn / a;
271 : :
272 [ - + ]: 1500 : if (Ma <= -1) { // Supersonic inflow
273 : : // For supersonic inflow, all the characteristics are from outside.
274 : : // Therefore, we calculate the ghost cell state using the primitive
275 : : // variables from outside.
276 : : tk::real rho(0.0);
277 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
278 [ - - ]: 0 : if (k == fmat)
279 : 0 : ur[volfracIdx(nmat,k)] = 1.0 -
280 : 0 : (static_cast< tk::real >(nmat-1))*alphamin;
281 : : else
282 : 0 : ur[volfracIdx(nmat,k)] = alphamin;
283 [ - - ]: 0 : auto rhok = mat_blk[k].compute< EOS::density >(fp, ft);
284 [ - - ]: 0 : ur[densityIdx(nmat,k)] = ur[volfracIdx(nmat,k)] * rhok;
285 : 0 : ur[energyIdx(nmat,k)] =
286 [ - - ]: 0 : mat_blk[k].compute< EOS::totalenergy >(ur[volfracIdx(nmat,k)]*rhok,
287 [ - - ]: 0 : fu[0], fu[1], fu[2], ur[volfracIdx(nmat,k)]*fp, ur[volfracIdx(nmat,k)]);
288 : :
289 : : // material pressures
290 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ur[volfracIdx(nmat, k)] * fp;
291 : :
292 : 0 : rho += ur[densityIdx(nmat,k)];
293 : : }
294 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
295 : 0 : ur[momentumIdx(nmat,i)] = rho * fu[i];
296 : 0 : ur[ncomp+velocityIdx(nmat, i)] = fu[i];
297 : : }
298 : :
299 [ + - ][ + - ]: 1500 : } else if (Ma > -1 && Ma < 0) { // Subsonic inflow
300 : : // For subsonic inflow, there is 1 outgoing characteristic and 4
301 : : // incoming characteristics. Therefore, we calculate the ghost cell state
302 : : // by taking pressure from the internal cell and other quantities from
303 : : // the outside.
304 : : tk::real rho(0.0);
305 [ + + ]: 4500 : for (std::size_t k=0; k<nmat; ++k) {
306 [ + + ]: 3000 : if (k == fmat)
307 : 1500 : ur[volfracIdx(nmat,k)] = 1.0 -
308 : 1500 : (static_cast< tk::real >(nmat-1))*alphamin;
309 : : else
310 : 1500 : ur[volfracIdx(nmat,k)] = alphamin;
311 [ + - ]: 3000 : auto p = ul[ncomp+pressureIdx(nmat,k)] / ul[volfracIdx(nmat,k)];
312 [ + - ]: 3000 : auto rhok = mat_blk[k].compute< EOS::density >(p, ft);
313 [ + - ]: 3000 : ur[densityIdx(nmat,k)] = ur[volfracIdx(nmat,k)] * rhok;
314 : 3000 : ur[energyIdx(nmat,k)] =
315 [ + - ]: 3000 : mat_blk[k].compute< EOS::totalenergy >(ur[volfracIdx(nmat,k)]*rhok,
316 [ + - ]: 3000 : fu[0], fu[1], fu[2], ur[volfracIdx(nmat,k)]*p, ur[volfracIdx(nmat,k)]);
317 : :
318 : : // material pressures
319 : 3000 : ur[ncomp+pressureIdx(nmat, k)] = ur[volfracIdx(nmat, k)] * p;
320 : :
321 : 3000 : rho += ur[densityIdx(nmat,k)];
322 : : }
323 [ + + ]: 6000 : for (std::size_t i=0; i<3; ++i) {
324 : 4500 : ur[momentumIdx(nmat,i)] = rho * fu[i];
325 : 4500 : ur[ncomp+velocityIdx(nmat, i)] = fu[i];
326 : : }
327 : :
328 [ - - ][ - - ]: 0 : } else if (Ma >= 0 && Ma < 1) { // Subsonic outflow
329 : : // For subsonic outflow, there is 1 incoming characteristic and 4
330 : : // outgoing characteristics. Therefore, we calculate the ghost cell state
331 : : // by taking pressure from the outside and other quantities from the
332 : : // internal cell.
333 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
334 : 0 : ur[energyIdx(nmat, k)] =
335 [ - - ]: 0 : mat_blk[k].compute< EOS::totalenergy >(
336 [ - - ]: 0 : ur[densityIdx(nmat, k)], v1l, v2l, v3l, ul[volfracIdx(nmat, k)]*fp,
337 [ - - ]: 0 : ul[volfracIdx(nmat, k)] );
338 : :
339 : : // material pressures
340 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ul[volfracIdx(nmat, k)] * fp;
341 : : }
342 : : }
343 : : // Otherwise, for supersonic outflow, all the characteristics are from
344 : : // internal cell. Therefore, we calculate the ghost cell state using the
345 : : // conservative variables from internal cell (which is what ur is
346 : : // initialized to).
347 : :
348 : : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
349 : : "boundary state vector" );
350 : :
351 [ + - ]: 3000 : return {{ std::move(ul), std::move(ur) }};
352 : : }
353 : :
354 : : //! \brief Boundary state function providing the left and right state of a
355 : : //! face at extrapolation boundaries
356 : : //! \param[in] ul Left (domain-internal) state
357 : : //! \return Left and right states for all scalar components in this PDE
358 : : //! system
359 : : //! \note The function signature must follow tk::StateFn. For multimat, the
360 : : //! left or right state is the vector of conserved quantities, followed by
361 : : //! the vector of primitive quantities appended to it.
362 : : static tk::StateFn::result_type
363 : 108976 : extrapolate( ncomp_t,
364 : : const std::vector< EOS >&,
365 : : const std::vector< tk::real >& ul,
366 : : tk::real, tk::real, tk::real, tk::real,
367 : : const std::array< tk::real, 3 >& )
368 : : {
369 : 108976 : return {{ ul, ul }};
370 : : }
371 : :
372 : : //! \brief Boundary state function providing the left and right state of a
373 : : //! face at no-slip wall boundaries
374 : : //! \param[in] ncomp Number of scalar components in this PDE system
375 : : //! \param[in] ul Left (domain-internal) state
376 : : //! \param[in] fn Unit face normal
377 : : //! \return Left and right states for all scalar components in this PDE
378 : : //! system
379 : : //! \note The function signature must follow tk::StateFn. For multimat, the
380 : : //! left or right state is the vector of conserved quantities, followed by
381 : : //! the vector of primitive quantities appended to it.
382 : : static tk::StateFn::result_type
383 : 0 : noslipwall( ncomp_t ncomp,
384 : : const std::vector< EOS >&,
385 : : const std::vector< tk::real >& ul,
386 : : tk::real, tk::real, tk::real, tk::real,
387 : : const std::array< tk::real, 3 >& fn )
388 : : {
389 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
390 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
391 : :
392 : 0 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
393 : :
394 : : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
395 : : "internal state vector" );
396 : :
397 : : tk::real rho(0.0);
398 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
399 : 0 : rho += ul[densityIdx(nmat, k)];
400 : :
401 : 0 : auto ur = ul;
402 : :
403 : : // Internal cell velocity components
404 : 0 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
405 : 0 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
406 : 0 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
407 : : // Ghost state velocity components
408 : 0 : auto v1r = -v1l;
409 : 0 : auto v2r = -v2l;
410 : 0 : auto v3r = -v3l;
411 : : // Boundary condition
412 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
413 : : {
414 [ - - ]: 0 : ur[volfracIdx(nmat, k)] = ul[volfracIdx(nmat, k)];
415 [ - - ]: 0 : ur[densityIdx(nmat, k)] = ul[densityIdx(nmat, k)];
416 : 0 : ur[energyIdx(nmat, k)] = ul[energyIdx(nmat, k)];
417 [ - - ]: 0 : if (solidx[k] > 0) {
418 : : // Internal inverse deformation tensor
419 : : std::array< std::array< tk::real, 3 >, 3 > g;
420 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
421 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
422 : 0 : g[i][j] = ul[deformIdx(nmat,solidx[k],i,j)];
423 : : // Internal Cauchy stress tensor
424 : : std::array< std::array< tk::real, 3 >, 3 > s;
425 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
426 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
427 : 0 : s[i][j] = ul[ncomp+stressIdx(nmat,solidx[k],stressCmp[i][j])];
428 : : // Make reflection matrix
429 : : std::array< std::array< tk::real, 3 >, 3 >
430 : 0 : reflectionMat{{{1,0,0}, {0,1,0}, {0,0,1}}};
431 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
432 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
433 : 0 : reflectionMat[i][j] -= 2*fn[i]*fn[j];
434 : : // Reflect g
435 [ - - ]: 0 : g = tk::reflectTensor(g, reflectionMat);
436 : : // Reflect s
437 [ - - ]: 0 : s = tk::reflectTensor(s, reflectionMat);
438 : : // Copy g and s into ur
439 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
440 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
441 : 0 : ur[deformIdx(nmat,solidx[k],i,j)] = g[i][j];
442 : 0 : ur[ncomp+stressIdx(nmat,solidx[k],stressCmp[i][j])] = s[i][j];
443 : : }
444 : : }
445 : : }
446 : 0 : ur[momentumIdx(nmat, 0)] = rho * v1r;
447 : 0 : ur[momentumIdx(nmat, 1)] = rho * v2r;
448 : 0 : ur[momentumIdx(nmat, 2)] = rho * v3r;
449 : :
450 : : // Internal cell primitive quantities using the separately reconstructed
451 : : // primitive quantities. This is used to get ghost state for primitive
452 : : // quantities
453 : :
454 : : // velocity
455 : 0 : ur[ncomp+velocityIdx(nmat, 0)] = v1r;
456 : 0 : ur[ncomp+velocityIdx(nmat, 1)] = v2r;
457 : 0 : ur[ncomp+velocityIdx(nmat, 2)] = v3r;
458 : : // material pressures
459 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
460 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ul[ncomp+pressureIdx(nmat, k)];
461 : :
462 : : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
463 : : "boundary state vector" );
464 : :
465 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
466 : : }
467 : :
468 : : //! \brief Boundary state function providing the left and right state of a
469 : : //! face at back pressure boundaries
470 : : //! \param[in] ncomp Number of scalar components in this PDE system
471 : : //! \param[in] ul Left (domain-internal) state
472 : : //! \return Left and right states for all scalar components in this PDE
473 : : //! system
474 : : //! \note The function signature must follow tk::StateFn
475 : : static tk::StateFn::result_type
476 : 0 : back_pressure( ncomp_t ncomp,
477 : : const std::vector< EOS >& mat_blk,
478 : : const std::vector< tk::real >& ul,
479 : : tk::real, tk::real, tk::real, tk::real,
480 : : const std::array< tk::real, 3 >& )
481 : : {
482 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
483 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
484 : :
485 : : // Back pressure
486 : 0 : auto fbp = g_inputdeck.get< tag::bc >()[0].get< tag::back_pressure >().get<
487 : 0 : tag::pressure >();
488 : :
489 : 0 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
490 : :
491 : : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
492 : : "internal state vector" );
493 : :
494 : 0 : auto ur = ul;
495 : :
496 : : // Internal cell velocity components
497 : : std::array< tk::real, 3 > vl;
498 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
499 : 0 : vl[i] = ul[ncomp+velocityIdx(nmat, i)];
500 : :
501 : : // The ghost cell state is calculated using the back pressure and other
502 : : // quantities from the internal cell
503 : : auto rhor_b(0.0);
504 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
505 [ - - ]: 0 : auto T_k = mat_blk[k].compute< inciter::EOS::temperature >(
506 : : ul[densityIdx(nmat, k)], vl[0], vl[1], vl[2], ul[energyIdx(nmat, k)],
507 [ - - ]: 0 : ul[volfracIdx(nmat, k)] );
508 [ - - ]: 0 : auto rhor_k = mat_blk[k].compute< inciter::EOS::density >( fbp, T_k );
509 [ - - ]: 0 : ur[densityIdx(nmat, k)] = ul[volfracIdx(nmat, k)] * rhor_k;
510 : 0 : ur[energyIdx(nmat, k)] =
511 : 0 : mat_blk[k].compute< EOS::totalenergy >(
512 : 0 : ul[volfracIdx(nmat, k)]*rhor_k, vl[0], vl[1], vl[2],
513 [ - - ][ - - ]: 0 : ul[volfracIdx(nmat, k)]*fbp, ul[volfracIdx(nmat, k)] );
514 : 0 : rhor_b += ur[densityIdx(nmat, k)];
515 : :
516 : : // material pressures
517 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ul[volfracIdx(nmat, k)] * fbp;
518 : : }
519 : :
520 : : // bulk momentum
521 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
522 : 0 : ur[momentumIdx(nmat,i)] = rhor_b * vl[i];
523 : :
524 : : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
525 : : "boundary state vector" );
526 : :
527 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
528 : : }
529 : :
530 : : //----------------------------------------------------------------------------
531 : : // Boundary Gradient functions
532 : : //----------------------------------------------------------------------------
533 : :
534 : : //! \brief Boundary gradient function copying the left gradient to the right
535 : : //! gradient at a face
536 : : //! \param[in] dul Left (domain-internal) state
537 : : //! \return Left and right states for all scalar components in this PDE
538 : : //! system
539 : : //! \note The function signature must follow tk::StateFn. For multimat, the
540 : : //! left or right state is the vector of gradients of primitive quantities.
541 : : static tk::StateFn::result_type
542 : 0 : noOpGrad( ncomp_t,
543 : : const std::vector< EOS >&,
544 : : const std::vector< tk::real >& dul,
545 : : tk::real, tk::real, tk::real, tk::real,
546 : : const std::array< tk::real, 3 >& )
547 : : {
548 : 0 : return {{ dul, dul }};
549 : : }
550 : :
551 : : //! \brief Boundary gradient function for the symmetry boundary condition
552 : : //! \param[in] ncomp Number of variables whos gradients are needed
553 : : //! \param[in] dul Left (domain-internal) gradients
554 : : //! \return Left and right states for all scalar components in this PDE
555 : : //! system
556 : : //! \note The function signature must follow tk::StateFn. For multimat, the
557 : : //! left or right state is the vector of gradients of primitive quantities.
558 : : static tk::StateFn::result_type
559 : 0 : symmetryGrad( ncomp_t ncomp,
560 : : const std::vector< EOS >&,
561 : : const std::vector< tk::real >& dul,
562 : : tk::real, tk::real, tk::real, tk::real,
563 : : const std::array< tk::real, 3 >& )
564 : : {
565 : : Assert(dul.size() == 3*ncomp, "Incorrect size of boundary gradient vector");
566 : :
567 : 0 : auto dur = dul;
568 : :
569 [ - - ]: 0 : for (std::size_t i=0; i<3*ncomp; ++i)
570 : 0 : dur[i] = -dul[i];
571 : :
572 [ - - ]: 0 : return {{ std::move(dul), std::move(dur) }};
573 : : }
574 : :
575 : : //! \brief Boundary gradient function for zero gradient cells
576 : : //! \param[in] ncomp Number of variables whos gradients are needed
577 : : //! \param[in] dul Left (domain-internal) gradients
578 : : //! \return Left and right states for all scalar components in this PDE
579 : : //! system
580 : : //! \note The function signature must follow tk::StateFn. For multimat, the
581 : : //! left or right state is the vector of gradients of primitive quantities.
582 : : static tk::StateFn::result_type
583 : 0 : zeroGrad( ncomp_t ncomp,
584 : : const std::vector< EOS >&,
585 : : const std::vector< tk::real >& dul,
586 : : tk::real, tk::real, tk::real, tk::real,
587 : : const std::array< tk::real, 3 >& )
588 : : {
589 : : Assert(dul.size() == 3*ncomp, "Incorrect size of boundary gradient vector");
590 : :
591 : 0 : std::vector< tk::real > dur(3*ncomp, 0.0);
592 : :
593 [ - - ]: 0 : return {{ std::move(dul), std::move(dur) }};
594 : : }
595 : :
596 : : } // inciter::
597 : :
598 : : #endif // BCFunctions_h
|
