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 : 1760240 : 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 : 1760240 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
39 : 1760240 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
40 : :
41 [ + - ]: 1760240 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
42 : :
43 [ - + ][ - - ]: 1760240 : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
44 : : "internal state vector" );
45 : :
46 : 1760240 : tk::real rho(0.0);
47 [ + + ]: 5904120 : for (std::size_t k=0; k<nmat; ++k)
48 : 4143880 : rho += ul[densityIdx(nmat, k)];
49 : :
50 [ + - ]: 3520480 : auto ur = ul;
51 : :
52 : : // Internal cell velocity components
53 : 1760240 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
54 : 1760240 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
55 : 1760240 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
56 : : // Normal component of velocity
57 : 1760240 : auto vnl = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
58 : : // Ghost state velocity components
59 : 1760240 : auto v1r = v1l - 2.0*vnl*fn[0];
60 : 1760240 : auto v2r = v2l - 2.0*vnl*fn[1];
61 : 1760240 : auto v3r = v3l - 2.0*vnl*fn[2];
62 : : // Boundary condition
63 [ + + ]: 5904120 : for (std::size_t k=0; k<nmat; ++k)
64 : : {
65 : 4143880 : ur[volfracIdx(nmat, k)] = ul[volfracIdx(nmat, k)];
66 : 4143880 : ur[densityIdx(nmat, k)] = ul[densityIdx(nmat, k)];
67 : 4143880 : ur[energyIdx(nmat, k)] = ul[energyIdx(nmat, k)];
68 [ - + ]: 4143880 : 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 : 1760240 : ur[momentumIdx(nmat, 0)] = rho * v1r;
98 : 1760240 : ur[momentumIdx(nmat, 1)] = rho * v2r;
99 : 1760240 : 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 : 1760240 : ur[ncomp+velocityIdx(nmat, 0)] = v1r;
107 : 1760240 : ur[ncomp+velocityIdx(nmat, 1)] = v2r;
108 : 1760240 : ur[ncomp+velocityIdx(nmat, 2)] = v3r;
109 : : // material pressures
110 [ + + ]: 5904120 : for (std::size_t k=0; k<nmat; ++k)
111 : 4143880 : ur[ncomp+pressureIdx(nmat, k)] = ul[ncomp+pressureIdx(nmat, k)];
112 : :
113 [ - + ][ - - ]: 1760240 : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
114 : : "boundary state vector" );
115 : :
116 [ + - ]: 3520480 : 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 : 0 : 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 [ - - ][ - - ]: 0 : 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 : 0 : ul[volfracIdx(nmat, k)], k ) );
168 : :
169 : : // Mach number
170 : 0 : auto Ma = vn / a;
171 : :
172 : 0 : tk::real alphamin = 1e-12;
173 : 0 : tk::real pk(0.0);
174 : 0 : 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)] = ur[volfracIdx(nmat,k)] *
193 [ - - ]: 0 : mat_blk[k].compute< EOS::totalenergy >(rhok, v1r, v2r, v3r, pk);
194 : :
195 : : // bulk density
196 : 0 : rho += ur[densityIdx(nmat,k)];
197 : : }
198 : :
199 : 0 : ur[momentumIdx(nmat, 0)] = rho * v1r;
200 : 0 : ur[momentumIdx(nmat, 1)] = rho * v2r;
201 : 0 : ur[momentumIdx(nmat, 2)] = rho * v3r;
202 : :
203 : : // velocity
204 : 0 : ur[ncomp+velocityIdx(nmat, 0)] = v1r;
205 : 0 : ur[ncomp+velocityIdx(nmat, 1)] = v2r;
206 : 0 : ur[ncomp+velocityIdx(nmat, 2)] = v3r;
207 : :
208 [ - - ][ - - ]: 0 : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
209 : : "boundary state vector" );
210 : :
211 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
212 : : }
213 : :
214 : : //! \brief Boundary state function providing the left and right state of a
215 : : //! face at farfield boundaries
216 : : //! \param[in] ncomp Number of scalar components in this PDE system
217 : : //! \param[in] ul Left (domain-internal) state
218 : : //! \param[in] fn Unit face normal
219 : : //! \return Left and right states for all scalar components in this PDE
220 : : //! system
221 : : //! \details The farfield boudary calculation, implemented here, is
222 : : //! based on the characteristic theory of hyperbolic systems.
223 : : //! \note The function signature must follow tk::StateFn
224 : : static tk::StateFn::result_type
225 : 0 : farfield( ncomp_t ncomp,
226 : : const std::vector< EOS >& mat_blk,
227 : : const std::vector< tk::real >& ul,
228 : : tk::real, tk::real, tk::real, tk::real,
229 : : const std::array< tk::real, 3 >& fn )
230 : : {
231 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
232 : 0 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
233 : :
234 : : // Farfield primitive quantities
235 : : auto fp =
236 : 0 : g_inputdeck.get< tag::bc >()[0].get< tag::pressure >();
237 : : auto ft =
238 : 0 : g_inputdeck.get< tag::bc >()[0].get< tag::temperature >();
239 : : auto fu =
240 [ - - ]: 0 : g_inputdeck.get< tag::bc >()[0].get< tag::velocity >();
241 : : auto fmat =
242 : 0 : g_inputdeck.get< tag::bc >()[0].get< tag::materialid >() - 1;
243 : :
244 [ - - ]: 0 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
245 : :
246 [ - - ][ - - ]: 0 : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
247 : : "internal state vector" );
248 : :
249 [ - - ]: 0 : auto ur = ul;
250 : :
251 : : // Internal cell velocity components
252 : 0 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
253 : 0 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
254 : 0 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
255 : :
256 : : // Normal velocity
257 : 0 : auto vn = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
258 : :
259 : : // Acoustic speed
260 : 0 : tk::real a(0.0);
261 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
262 [ - - ]: 0 : if (ul[volfracIdx(nmat, k)] > 1.0e-04)
263 [ - - ]: 0 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
264 : 0 : ul[densityIdx(nmat, k)], ul[ncomp+pressureIdx(nmat, k)],
265 : 0 : ul[volfracIdx(nmat, k)], k ) );
266 : :
267 : : // Mach number
268 : 0 : auto Ma = vn / a;
269 : :
270 : 0 : tk::real alphamin = 1e-12;
271 : :
272 [ - - ]: 0 : 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 : 0 : 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)] = ur[volfracIdx(nmat,k)] *
286 [ - - ]: 0 : mat_blk[k].compute< EOS::totalenergy >(rhok, fu[0], fu[1], fu[2], fp);
287 : :
288 : : // material pressures
289 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ur[volfracIdx(nmat, k)] * fp;
290 : :
291 : 0 : rho += ur[densityIdx(nmat,k)];
292 : : }
293 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
294 : 0 : ur[momentumIdx(nmat,i)] = rho * fu[i];
295 : 0 : ur[ncomp+velocityIdx(nmat, i)] = fu[i];
296 : : }
297 : :
298 [ - - ][ - - ]: 0 : } else if (Ma > -1 && Ma < 0) { // Subsonic inflow
299 : : // For subsonic inflow, there is 1 outgoing characteristic and 4
300 : : // incoming characteristics. Therefore, we calculate the ghost cell state
301 : : // by taking pressure from the internal cell and other quantities from
302 : : // the outside.
303 : 0 : tk::real rho(0.0);
304 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
305 [ - - ]: 0 : if (k == fmat)
306 : 0 : ur[volfracIdx(nmat,k)] = 1.0 -
307 : 0 : (static_cast< tk::real >(nmat-1))*alphamin;
308 : : else
309 : 0 : ur[volfracIdx(nmat,k)] = alphamin;
310 : 0 : auto p = ul[ncomp+pressureIdx(nmat,k)] / ul[volfracIdx(nmat,k)];
311 [ - - ]: 0 : auto rhok = mat_blk[k].compute< EOS::density >(p, ft);
312 : 0 : ur[densityIdx(nmat,k)] = ur[volfracIdx(nmat,k)] * rhok;
313 : 0 : ur[energyIdx(nmat,k)] = ur[volfracIdx(nmat,k)] *
314 [ - - ]: 0 : mat_blk[k].compute< EOS::totalenergy >(rhok, fu[0], fu[1], fu[2], p);
315 : :
316 : : // material pressures
317 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ur[volfracIdx(nmat, k)] * p;
318 : :
319 : 0 : rho += ur[densityIdx(nmat,k)];
320 : : }
321 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
322 : 0 : ur[momentumIdx(nmat,i)] = rho * fu[i];
323 : 0 : ur[ncomp+velocityIdx(nmat, i)] = fu[i];
324 : 0 : }
325 : :
326 [ - - ][ - - ]: 0 : } else if (Ma >= 0 && Ma < 1) { // Subsonic outflow
327 : : // For subsonic outflow, there is 1 incoming characteristic and 4
328 : : // outgoing characteristics. Therefore, we calculate the ghost cell state
329 : : // by taking pressure from the outside and other quantities from the
330 : : // internal cell.
331 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
332 : 0 : ur[energyIdx(nmat, k)] = ul[volfracIdx(nmat, k)] *
333 : 0 : mat_blk[k].compute< EOS::totalenergy >(
334 [ - - ]: 0 : ur[densityIdx(nmat, k)]/ul[volfracIdx(nmat, k)], v1l, v2l, v3l, fp );
335 : :
336 : : // material pressures
337 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ul[volfracIdx(nmat, k)] * fp;
338 : : }
339 : : }
340 : : // Otherwise, for supersonic outflow, all the characteristics are from
341 : : // internal cell. Therefore, we calculate the ghost cell state using the
342 : : // conservative variables from internal cell (which is what ur is
343 : : // initialized to).
344 : :
345 [ - - ][ - - ]: 0 : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
346 : : "boundary state vector" );
347 : :
348 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
349 : : }
350 : :
351 : : //! \brief Boundary state function providing the left and right state of a
352 : : //! face at extrapolation boundaries
353 : : //! \param[in] ul Left (domain-internal) state
354 : : //! \return Left and right states for all scalar components in this PDE
355 : : //! system
356 : : //! \note The function signature must follow tk::StateFn. For multimat, the
357 : : //! left or right state is the vector of conserved quantities, followed by
358 : : //! the vector of primitive quantities appended to it.
359 : : static tk::StateFn::result_type
360 : 106676 : extrapolate( ncomp_t,
361 : : const std::vector< EOS >&,
362 : : const std::vector< tk::real >& ul,
363 : : tk::real, tk::real, tk::real, tk::real,
364 : : const std::array< tk::real, 3 >& )
365 : : {
366 : 106676 : return {{ ul, ul }};
367 : : }
368 : :
369 : : //! \brief Boundary state function providing the left and right state of a
370 : : //! face at no-slip wall boundaries
371 : : //! \param[in] ncomp Number of scalar components in this PDE system
372 : : //! \param[in] ul Left (domain-internal) state
373 : : //! \param[in] fn Unit face normal
374 : : //! \return Left and right states for all scalar components in this PDE
375 : : //! system
376 : : //! \note The function signature must follow tk::StateFn. For multimat, the
377 : : //! left or right state is the vector of conserved quantities, followed by
378 : : //! the vector of primitive quantities appended to it.
379 : : static tk::StateFn::result_type
380 : 0 : noslipwall( ncomp_t ncomp,
381 : : const std::vector< EOS >&,
382 : : const std::vector< tk::real >& ul,
383 : : tk::real, tk::real, tk::real, tk::real,
384 : : const std::array< tk::real, 3 >& fn )
385 : : {
386 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
387 : 0 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
388 : :
389 [ - - ]: 0 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
390 : :
391 [ - - ][ - - ]: 0 : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
392 : : "internal state vector" );
393 : :
394 : 0 : tk::real rho(0.0);
395 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
396 : 0 : rho += ul[densityIdx(nmat, k)];
397 : :
398 [ - - ]: 0 : auto ur = ul;
399 : :
400 : : // Internal cell velocity components
401 : 0 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
402 : 0 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
403 : 0 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
404 : : // Ghost state velocity components
405 : 0 : auto v1r = -v1l;
406 : 0 : auto v2r = -v2l;
407 : 0 : auto v3r = -v3l;
408 : : // Boundary condition
409 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
410 : : {
411 : 0 : ur[volfracIdx(nmat, k)] = ul[volfracIdx(nmat, k)];
412 : 0 : ur[densityIdx(nmat, k)] = ul[densityIdx(nmat, k)];
413 : 0 : ur[energyIdx(nmat, k)] = ul[energyIdx(nmat, k)];
414 [ - - ]: 0 : if (solidx[k] > 0) {
415 : : // Internal inverse deformation tensor
416 : : std::array< std::array< tk::real, 3 >, 3 > g;
417 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
418 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
419 : 0 : g[i][j] = ul[deformIdx(nmat,solidx[k],i,j)];
420 : : // Internal Cauchy stress tensor
421 : : std::array< std::array< tk::real, 3 >, 3 > s;
422 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
423 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
424 : 0 : s[i][j] = ul[ncomp+stressIdx(nmat,solidx[k],stressCmp[i][j])];
425 : : // Make reflection matrix
426 : : std::array< std::array< tk::real, 3 >, 3 >
427 : 0 : reflectionMat{{{1,0,0}, {0,1,0}, {0,0,1}}};
428 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
429 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
430 : 0 : reflectionMat[i][j] -= 2*fn[i]*fn[j];
431 : : // Reflect g
432 [ - - ]: 0 : g = tk::reflectTensor(g, reflectionMat);
433 : : // Reflect s
434 [ - - ]: 0 : s = tk::reflectTensor(s, reflectionMat);
435 : : // Copy g and s into ur
436 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
437 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
438 : 0 : ur[deformIdx(nmat,solidx[k],i,j)] = g[i][j];
439 : 0 : ur[ncomp+stressIdx(nmat,solidx[k],stressCmp[i][j])] = s[i][j];
440 : : }
441 : : }
442 : : }
443 : 0 : ur[momentumIdx(nmat, 0)] = rho * v1r;
444 : 0 : ur[momentumIdx(nmat, 1)] = rho * v2r;
445 : 0 : ur[momentumIdx(nmat, 2)] = rho * v3r;
446 : :
447 : : // Internal cell primitive quantities using the separately reconstructed
448 : : // primitive quantities. This is used to get ghost state for primitive
449 : : // quantities
450 : :
451 : : // velocity
452 : 0 : ur[ncomp+velocityIdx(nmat, 0)] = v1r;
453 : 0 : ur[ncomp+velocityIdx(nmat, 1)] = v2r;
454 : 0 : ur[ncomp+velocityIdx(nmat, 2)] = v3r;
455 : : // material pressures
456 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
457 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ul[ncomp+pressureIdx(nmat, k)];
458 : :
459 [ - - ][ - - ]: 0 : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
460 : : "boundary state vector" );
461 : :
462 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
463 : : }
464 : :
465 : : //----------------------------------------------------------------------------
466 : : // Boundary Gradient functions
467 : : //----------------------------------------------------------------------------
468 : :
469 : : //! \brief Boundary gradient function copying the left gradient to the right
470 : : //! gradient at a face
471 : : //! \param[in] dul 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. For multimat, the
475 : : //! left or right state is the vector of gradients of primitive quantities.
476 : : static tk::StateFn::result_type
477 : 0 : noOpGrad( ncomp_t,
478 : : const std::vector< EOS >&,
479 : : const std::vector< tk::real >& dul,
480 : : tk::real, tk::real, tk::real, tk::real,
481 : : const std::array< tk::real, 3 >& )
482 : : {
483 : 0 : return {{ dul, dul }};
484 : : }
485 : :
486 : : //! \brief Boundary gradient function for the symmetry boundary condition
487 : : //! \param[in] ncomp Number of variables whos gradients are needed
488 : : //! \param[in] dul Left (domain-internal) gradients
489 : : //! \return Left and right states for all scalar components in this PDE
490 : : //! system
491 : : //! \note The function signature must follow tk::StateFn. For multimat, the
492 : : //! left or right state is the vector of gradients of primitive quantities.
493 : : static tk::StateFn::result_type
494 : 0 : symmetryGrad( ncomp_t ncomp,
495 : : const std::vector< EOS >&,
496 : : const std::vector< tk::real >& dul,
497 : : tk::real, tk::real, tk::real, tk::real,
498 : : const std::array< tk::real, 3 >& )
499 : : {
500 [ - - ][ - - ]: 0 : Assert(dul.size() == 3*ncomp, "Incorrect size of boundary gradient vector");
[ - - ][ - - ]
501 : :
502 [ - - ]: 0 : auto dur = dul;
503 : :
504 [ - - ]: 0 : for (std::size_t i=0; i<3*ncomp; ++i)
505 : 0 : dur[i] = -dul[i];
506 : :
507 [ - - ]: 0 : return {{ std::move(dul), std::move(dur) }};
508 : : }
509 : :
510 : : //! \brief Boundary gradient function for zero gradient cells
511 : : //! \param[in] ncomp Number of variables whos gradients are needed
512 : : //! \param[in] dul Left (domain-internal) gradients
513 : : //! \return Left and right states for all scalar components in this PDE
514 : : //! system
515 : : //! \note The function signature must follow tk::StateFn. For multimat, the
516 : : //! left or right state is the vector of gradients of primitive quantities.
517 : : static tk::StateFn::result_type
518 : 0 : zeroGrad( ncomp_t ncomp,
519 : : const std::vector< EOS >&,
520 : : const std::vector< tk::real >& dul,
521 : : tk::real, tk::real, tk::real, tk::real,
522 : : const std::array< tk::real, 3 >& )
523 : : {
524 [ - - ][ - - ]: 0 : Assert(dul.size() == 3*ncomp, "Incorrect size of boundary gradient vector");
[ - - ][ - - ]
525 : :
526 [ - - ]: 0 : std::vector< tk::real > dur(3*ncomp, 0.0);
527 : :
528 [ - - ]: 0 : return {{ std::move(dul), std::move(dur) }};
529 : : }
530 : :
531 : : } // inciter::
532 : :
533 : : #endif // BCFunctions_h
|
