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 : 1818040 : 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 : 1818040 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
39 : 1818040 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
40 : :
41 [ + - ]: 1818040 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
42 : :
43 [ - + ][ - - ]: 1818040 : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
44 : : "internal state vector" );
45 : :
46 : 1818040 : tk::real rho(0.0);
47 [ + + ]: 6077520 : for (std::size_t k=0; k<nmat; ++k)
48 : 4259480 : rho += ul[densityIdx(nmat, k)];
49 : :
50 [ + - ]: 3636080 : auto ur = ul;
51 : :
52 : : // Internal cell velocity components
53 : 1818040 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
54 : 1818040 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
55 : 1818040 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
56 : : // Normal component of velocity
57 : 1818040 : auto vnl = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
58 : : // Ghost state velocity components
59 : 1818040 : auto v1r = v1l - 2.0*vnl*fn[0];
60 : 1818040 : auto v2r = v2l - 2.0*vnl*fn[1];
61 : 1818040 : auto v3r = v3l - 2.0*vnl*fn[2];
62 : : // Boundary condition
63 [ + + ]: 6077520 : for (std::size_t k=0; k<nmat; ++k)
64 : : {
65 : 4259480 : ur[volfracIdx(nmat, k)] = ul[volfracIdx(nmat, k)];
66 : 4259480 : ur[densityIdx(nmat, k)] = ul[densityIdx(nmat, k)];
67 : 4259480 : ur[energyIdx(nmat, k)] = ul[energyIdx(nmat, k)];
68 [ - + ]: 4259480 : 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 : 1818040 : ur[momentumIdx(nmat, 0)] = rho * v1r;
98 : 1818040 : ur[momentumIdx(nmat, 1)] = rho * v2r;
99 : 1818040 : 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 : 1818040 : ur[ncomp+velocityIdx(nmat, 0)] = v1r;
107 : 1818040 : ur[ncomp+velocityIdx(nmat, 1)] = v2r;
108 : 1818040 : ur[ncomp+velocityIdx(nmat, 2)] = v3r;
109 : : // material pressures
110 [ + + ]: 6077520 : for (std::size_t k=0; k<nmat; ++k)
111 : 4259480 : ur[ncomp+pressureIdx(nmat, k)] = ul[ncomp+pressureIdx(nmat, k)];
112 : :
113 [ - + ][ - - ]: 1818040 : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
114 : : "boundary state vector" );
115 : :
116 [ + - ]: 3636080 : 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 farfield 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 farfield boudary calculation, implemented here, is
127 : : //! based on the characteristic theory of hyperbolic systems.
128 : : //! \note The function signature must follow tk::StateFn
129 : : static tk::StateFn::result_type
130 : 0 : farfield( 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 : :
139 : : // Farfield primitive quantities
140 : : auto fp =
141 : 0 : g_inputdeck.get< tag::bc >()[0].get< tag::pressure >();
142 : : auto ft =
143 : 0 : g_inputdeck.get< tag::bc >()[0].get< tag::temperature >();
144 : : auto fu =
145 [ - - ]: 0 : g_inputdeck.get< tag::bc >()[0].get< tag::velocity >();
146 : : auto fmat =
147 : 0 : g_inputdeck.get< tag::bc >()[0].get< tag::materialid >() - 1;
148 : :
149 [ - - ]: 0 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
150 : :
151 [ - - ][ - - ]: 0 : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
152 : : "internal state vector" );
153 : :
154 [ - - ]: 0 : auto ur = ul;
155 : :
156 : : // Internal cell velocity components
157 : 0 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
158 : 0 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
159 : 0 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
160 : :
161 : : // Normal velocity
162 : 0 : auto vn = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
163 : :
164 : : // Acoustic speed
165 : 0 : tk::real a(0.0);
166 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
167 [ - - ]: 0 : if (ul[volfracIdx(nmat, k)] > 1.0e-04)
168 [ - - ]: 0 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
169 : 0 : ul[densityIdx(nmat, k)], ul[ncomp+pressureIdx(nmat, k)],
170 : 0 : ul[volfracIdx(nmat, k)], k ) );
171 : :
172 : : // Mach number
173 : 0 : auto Ma = vn / a;
174 : :
175 : 0 : tk::real alphamin = 1e-12;
176 : :
177 [ - - ]: 0 : if (Ma <= -1) { // Supersonic inflow
178 : : // For supersonic inflow, all the characteristics are from outside.
179 : : // Therefore, we calculate the ghost cell state using the primitive
180 : : // variables from outside.
181 : 0 : tk::real rho(0.0);
182 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
183 [ - - ]: 0 : if (k == fmat)
184 : 0 : ur[volfracIdx(nmat,k)] = 1.0 -
185 : 0 : (static_cast< tk::real >(nmat-1))*alphamin;
186 : : else
187 : 0 : ur[volfracIdx(nmat,k)] = alphamin;
188 [ - - ]: 0 : auto rhok = mat_blk[k].compute< EOS::density >(fp, ft);
189 : 0 : ur[densityIdx(nmat,k)] = ur[volfracIdx(nmat,k)] * rhok;
190 : 0 : ur[energyIdx(nmat,k)] = ur[volfracIdx(nmat,k)] *
191 [ - - ]: 0 : mat_blk[k].compute< EOS::totalenergy >(rhok, fu[0], fu[1], fu[2], fp);
192 : :
193 : : // material pressures
194 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ul[volfracIdx(nmat, k)] * fp;
195 : :
196 : 0 : rho += ur[densityIdx(nmat,k)];
197 : : }
198 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
199 : 0 : ur[momentumIdx(nmat,i)] = rho * fu[i];
200 : 0 : ur[ncomp+velocityIdx(nmat, i)] = fu[i];
201 : : }
202 : :
203 [ - - ][ - - ]: 0 : } else if (Ma > -1 && Ma < 0) { // Subsonic inflow
204 : : // For subsonic inflow, there is 1 outgoing characteristic and 4
205 : : // incoming characteristics. Therefore, we calculate the ghost cell state
206 : : // by taking pressure from the internal cell and other quantities from
207 : : // the outside.
208 : 0 : tk::real rho(0.0);
209 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
210 [ - - ]: 0 : if (k == fmat)
211 : 0 : ur[volfracIdx(nmat,k)] = 1.0 -
212 : 0 : (static_cast< tk::real >(nmat-1))*alphamin;
213 : : else
214 : 0 : ur[volfracIdx(nmat,k)] = alphamin;
215 : 0 : auto p = ul[ncomp+pressureIdx(nmat,k)] / ul[volfracIdx(nmat,k)];
216 [ - - ]: 0 : auto rhok = mat_blk[k].compute< EOS::density >(p, ft);
217 : 0 : ur[densityIdx(nmat,k)] = ur[volfracIdx(nmat,k)] * rhok;
218 : 0 : ur[energyIdx(nmat,k)] = ur[volfracIdx(nmat,k)] *
219 [ - - ]: 0 : mat_blk[k].compute< EOS::totalenergy >(rhok, fu[0], fu[1], fu[2], p);
220 : :
221 : : // material pressures
222 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ul[volfracIdx(nmat, k)] * p;
223 : :
224 : 0 : rho += ur[densityIdx(nmat,k)];
225 : : }
226 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
227 : 0 : ur[momentumIdx(nmat,i)] = rho * fu[i];
228 : 0 : ur[ncomp+velocityIdx(nmat, i)] = fu[i];
229 : 0 : }
230 : :
231 [ - - ][ - - ]: 0 : } else if (Ma >= 0 && Ma < 1) { // Subsonic outflow
232 : : // For subsonic outflow, there is 1 incoming characteristic and 4
233 : : // outgoing characteristics. Therefore, we calculate the ghost cell state
234 : : // by taking pressure from the outside and other quantities from the
235 : : // internal cell.
236 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
237 : 0 : ur[energyIdx(nmat, k)] = ul[volfracIdx(nmat, k)] *
238 : 0 : mat_blk[k].compute< EOS::totalenergy >(
239 [ - - ]: 0 : ur[densityIdx(nmat, k)]/ul[volfracIdx(nmat, k)], v1l, v2l, v3l, fp );
240 : :
241 : : // material pressures
242 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ul[volfracIdx(nmat, k)] * fp;
243 : : }
244 : : }
245 : : // Otherwise, for supersonic outflow, all the characteristics are from
246 : : // internal cell. Therefore, we calculate the ghost cell state using the
247 : : // conservative variables from internal cell (which is what ur is
248 : : // initialized to).
249 : :
250 [ - - ][ - - ]: 0 : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
251 : : "boundary state vector" );
252 : :
253 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
254 : : }
255 : :
256 : : //! \brief Boundary state function providing the left and right state of a
257 : : //! face at extrapolation boundaries
258 : : //! \param[in] ul Left (domain-internal) state
259 : : //! \return Left and right states for all scalar components in this PDE
260 : : //! system
261 : : //! \note The function signature must follow tk::StateFn. For multimat, the
262 : : //! left or right state is the vector of conserved quantities, followed by
263 : : //! the vector of primitive quantities appended to it.
264 : : static tk::StateFn::result_type
265 : 106676 : extrapolate( ncomp_t,
266 : : const std::vector< EOS >&,
267 : : const std::vector< tk::real >& ul,
268 : : tk::real, tk::real, tk::real, tk::real,
269 : : const std::array< tk::real, 3 >& )
270 : : {
271 : 106676 : return {{ ul, ul }};
272 : : }
273 : :
274 : : //! \brief Boundary state function providing the left and right state of a
275 : : //! face at no-slip wall boundaries
276 : : //! \param[in] ncomp Number of scalar components in this PDE system
277 : : //! \param[in] ul Left (domain-internal) state
278 : : //! \param[in] fn Unit face normal
279 : : //! \return Left and right states for all scalar components in this PDE
280 : : //! system
281 : : //! \note The function signature must follow tk::StateFn. For multimat, the
282 : : //! left or right state is the vector of conserved quantities, followed by
283 : : //! the vector of primitive quantities appended to it.
284 : : static tk::StateFn::result_type
285 : 0 : noslipwall( ncomp_t ncomp,
286 : : const std::vector< EOS >&,
287 : : const std::vector< tk::real >& ul,
288 : : tk::real, tk::real, tk::real, tk::real,
289 : : const std::array< tk::real, 3 >& fn )
290 : : {
291 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
292 : 0 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
293 : :
294 [ - - ]: 0 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
295 : :
296 [ - - ][ - - ]: 0 : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
297 : : "internal state vector" );
298 : :
299 : 0 : tk::real rho(0.0);
300 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
301 : 0 : rho += ul[densityIdx(nmat, k)];
302 : :
303 [ - - ]: 0 : auto ur = ul;
304 : :
305 : : // Internal cell velocity components
306 : 0 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
307 : 0 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
308 : 0 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
309 : : // Ghost state velocity components
310 : 0 : auto v1r = -v1l;
311 : 0 : auto v2r = -v2l;
312 : 0 : auto v3r = -v3l;
313 : : // Boundary condition
314 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
315 : : {
316 : 0 : ur[volfracIdx(nmat, k)] = ul[volfracIdx(nmat, k)];
317 : 0 : ur[densityIdx(nmat, k)] = ul[densityIdx(nmat, k)];
318 : 0 : ur[energyIdx(nmat, k)] = ul[energyIdx(nmat, k)];
319 [ - - ]: 0 : if (solidx[k] > 0) {
320 : : // Internal inverse deformation tensor
321 : : std::array< std::array< tk::real, 3 >, 3 > g;
322 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
323 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
324 : 0 : g[i][j] = ul[deformIdx(nmat,solidx[k],i,j)];
325 : : // Internal Cauchy stress tensor
326 : : std::array< std::array< tk::real, 3 >, 3 > s;
327 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
328 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
329 : 0 : s[i][j] = ul[ncomp+stressIdx(nmat,solidx[k],stressCmp[i][j])];
330 : : // Make reflection matrix
331 : : std::array< std::array< tk::real, 3 >, 3 >
332 : 0 : reflectionMat{{{1,0,0}, {0,1,0}, {0,0,1}}};
333 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
334 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
335 : 0 : reflectionMat[i][j] -= 2*fn[i]*fn[j];
336 : : // Reflect g
337 [ - - ]: 0 : g = tk::reflectTensor(g, reflectionMat);
338 : : // Reflect s
339 [ - - ]: 0 : s = tk::reflectTensor(s, reflectionMat);
340 : : // Copy g and s into ur
341 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
342 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
343 : 0 : ur[deformIdx(nmat,solidx[k],i,j)] = g[i][j];
344 : 0 : ur[ncomp+stressIdx(nmat,solidx[k],stressCmp[i][j])] = s[i][j];
345 : : }
346 : : }
347 : : }
348 : 0 : ur[momentumIdx(nmat, 0)] = rho * v1r;
349 : 0 : ur[momentumIdx(nmat, 1)] = rho * v2r;
350 : 0 : ur[momentumIdx(nmat, 2)] = rho * v3r;
351 : :
352 : : // Internal cell primitive quantities using the separately reconstructed
353 : : // primitive quantities. This is used to get ghost state for primitive
354 : : // quantities
355 : :
356 : : // velocity
357 : 0 : ur[ncomp+velocityIdx(nmat, 0)] = v1r;
358 : 0 : ur[ncomp+velocityIdx(nmat, 1)] = v2r;
359 : 0 : ur[ncomp+velocityIdx(nmat, 2)] = v3r;
360 : : // material pressures
361 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
362 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ul[ncomp+pressureIdx(nmat, k)];
363 : :
364 [ - - ][ - - ]: 0 : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
365 : : "boundary state vector" );
366 : :
367 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
368 : : }
369 : :
370 : : //----------------------------------------------------------------------------
371 : : // Boundary Gradient functions
372 : : //----------------------------------------------------------------------------
373 : :
374 : : //! \brief Boundary gradient function copying the left gradient to the right
375 : : //! gradient at a face
376 : : //! \param[in] dul Left (domain-internal) state
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 gradients of primitive quantities.
381 : : static tk::StateFn::result_type
382 : 0 : noOpGrad( ncomp_t,
383 : : const std::vector< EOS >&,
384 : : const std::vector< tk::real >& dul,
385 : : tk::real, tk::real, tk::real, tk::real,
386 : : const std::array< tk::real, 3 >& )
387 : : {
388 : 0 : return {{ dul, dul }};
389 : : }
390 : :
391 : : //! \brief Boundary gradient function for the symmetry boundary condition
392 : : //! \param[in] ncomp Number of variables whos gradients are needed
393 : : //! \param[in] dul Left (domain-internal) gradients
394 : : //! \return Left and right states for all scalar components in this PDE
395 : : //! system
396 : : //! \note The function signature must follow tk::StateFn. For multimat, the
397 : : //! left or right state is the vector of gradients of primitive quantities.
398 : : static tk::StateFn::result_type
399 : 0 : symmetryGrad( ncomp_t ncomp,
400 : : const std::vector< EOS >&,
401 : : const std::vector< tk::real >& dul,
402 : : tk::real, tk::real, tk::real, tk::real,
403 : : const std::array< tk::real, 3 >& )
404 : : {
405 [ - - ][ - - ]: 0 : Assert(dul.size() == 3*ncomp, "Incorrect size of boundary gradient vector");
[ - - ][ - - ]
406 : :
407 [ - - ]: 0 : auto dur = dul;
408 : :
409 [ - - ]: 0 : for (std::size_t i=0; i<3*ncomp; ++i)
410 : 0 : dur[i] = -dul[i];
411 : :
412 [ - - ]: 0 : return {{ std::move(dul), std::move(dur) }};
413 : : }
414 : :
415 : : } // inciter::
416 : :
417 : : #endif // BCFunctions_h
|
