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 : 1876080 : 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 : 1876080 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
39 : 1876080 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
40 : :
41 [ + - ]: 1876080 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
42 : :
43 [ - + ][ - - ]: 1876080 : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
44 : : "internal state vector" );
45 : :
46 : 1876080 : tk::real rho(0.0);
47 [ + + ]: 6251640 : for (std::size_t k=0; k<nmat; ++k)
48 : 4375560 : rho += ul[densityIdx(nmat, k)];
49 : :
50 [ + - ]: 3752160 : auto ur = ul;
51 : :
52 : : // Internal cell velocity components
53 : 1876080 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
54 : 1876080 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
55 : 1876080 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
56 : : // Normal component of velocity
57 : 1876080 : auto vnl = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
58 : : // Ghost state velocity components
59 : 1876080 : auto v1r = v1l - 2.0*vnl*fn[0];
60 : 1876080 : auto v2r = v2l - 2.0*vnl*fn[1];
61 : 1876080 : auto v3r = v3l - 2.0*vnl*fn[2];
62 : : // Boundary condition
63 [ + + ]: 6251640 : for (std::size_t k=0; k<nmat; ++k)
64 : : {
65 : 4375560 : ur[volfracIdx(nmat, k)] = ul[volfracIdx(nmat, k)];
66 : 4375560 : ur[densityIdx(nmat, k)] = ul[densityIdx(nmat, k)];
67 : 4375560 : ur[energyIdx(nmat, k)] = ul[energyIdx(nmat, k)];
68 [ + + ]: 4375560 : if (solidx[k] > 0) {
69 : : // Internal inverse deformation tensor
70 : : std::array< std::array< tk::real, 3 >, 3 > g;
71 [ + + ]: 216000 : for (std::size_t i=0; i<3; ++i)
72 [ + + ]: 648000 : for (std::size_t j=0; j<3; ++j)
73 : 486000 : g[i][j] = ul[deformIdx(nmat,solidx[k],i,j)];
74 : : // Make reflection matrix
75 : : std::array< std::array< tk::real, 3 >, 3 >
76 : 54000 : reflectionMat{{{1,0,0}, {0,1,0}, {0,0,1}}};
77 [ + + ]: 216000 : for (std::size_t i=0; i<3; ++i)
78 [ + + ]: 648000 : for (std::size_t j=0; j<3; ++j)
79 : 486000 : reflectionMat[i][j] -= 2*fn[i]*fn[j];
80 : : // Reflect g
81 [ + - ]: 54000 : g = tk::reflectTensor(g, reflectionMat);
82 : : // Copy g into ur
83 [ + + ]: 216000 : for (std::size_t i=0; i<3; ++i)
84 [ + + ]: 648000 : for (std::size_t j=0; j<3; ++j)
85 : 486000 : ur[deformIdx(nmat,solidx[k],i,j)] = g[i][j];
86 : : }
87 : : }
88 : 1876080 : ur[momentumIdx(nmat, 0)] = rho * v1r;
89 : 1876080 : ur[momentumIdx(nmat, 1)] = rho * v2r;
90 : 1876080 : ur[momentumIdx(nmat, 2)] = rho * v3r;
91 : :
92 : : // Internal cell primitive quantities using the separately reconstructed
93 : : // primitive quantities. This is used to get ghost state for primitive
94 : : // quantities
95 : :
96 : : // velocity
97 : 1876080 : ur[ncomp+velocityIdx(nmat, 0)] = v1r;
98 : 1876080 : ur[ncomp+velocityIdx(nmat, 1)] = v2r;
99 : 1876080 : ur[ncomp+velocityIdx(nmat, 2)] = v3r;
100 : : // material pressures
101 [ + + ]: 6251640 : for (std::size_t k=0; k<nmat; ++k)
102 : 4375560 : ur[ncomp+pressureIdx(nmat, k)] = ul[ncomp+pressureIdx(nmat, k)];
103 : :
104 [ - + ][ - - ]: 1876080 : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
105 : : "boundary state vector" );
106 : :
107 [ + - ]: 3752160 : return {{ std::move(ul), std::move(ur) }};
108 : : }
109 : :
110 : : //! \brief Boundary state function providing the left and right state of a
111 : : //! face at farfield outlet boundaries
112 : : //! \param[in] ncomp Number of scalar components in this PDE system
113 : : //! \param[in] ul Left (domain-internal) state
114 : : //! \param[in] fn Unit face normal
115 : : //! \return Left and right states for all scalar components in this PDE
116 : : //! system
117 : : //! \details The farfield outlet boudary calculation, implemented here, is
118 : : //! based on the characteristic theory of hyperbolic systems. For subsonic
119 : : //! outlet flow, there is 1 incoming characteristic per material.
120 : : //! Therefore, we calculate the ghost cell state by taking material
121 : : //! pressure from the outside and other quantities from the internal cell.
122 : : //! For supersonic outlet flow, all the characteristics are from internal
123 : : //! cell and we obtain the ghost cell state from the internal cell.
124 : : //! \note The function signature must follow tk::StateFn
125 : : static tk::StateFn::result_type
126 : 0 : farfieldOutlet( ncomp_t ncomp,
127 : : const std::vector< EOS >& mat_blk,
128 : : const std::vector< tk::real >& ul,
129 : : tk::real, tk::real, tk::real, tk::real,
130 : : const std::array< tk::real, 3 >& fn )
131 : : {
132 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
133 : 0 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
134 : :
135 : : auto fp =
136 : 0 : g_inputdeck.get< tag::bc >()[0].get< tag::pressure >();
137 : :
138 [ - - ]: 0 : [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
139 : :
140 [ - - ][ - - ]: 0 : Assert( ul.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
141 : : "internal state vector" );
142 : :
143 [ - - ]: 0 : auto ur = ul;
144 : :
145 : : // Internal cell velocity components
146 : 0 : auto v1l = ul[ncomp+velocityIdx(nmat, 0)];
147 : 0 : auto v2l = ul[ncomp+velocityIdx(nmat, 1)];
148 : 0 : auto v3l = ul[ncomp+velocityIdx(nmat, 2)];
149 : :
150 : : // Normal velocity
151 : 0 : auto vn = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
152 : :
153 : : // Acoustic speed
154 : 0 : tk::real a(0.0);
155 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
156 [ - - ]: 0 : if (ul[volfracIdx(nmat, k)] > 1.0e-04)
157 [ - - ]: 0 : a = std::max( a, mat_blk[k].compute< EOS::soundspeed >(
158 : 0 : ul[densityIdx(nmat, k)], ul[ncomp+pressureIdx(nmat, k)],
159 : 0 : ul[volfracIdx(nmat, k)], k ) );
160 : :
161 : : // Mach number
162 : 0 : auto Ma = vn / a;
163 : :
164 [ - - ][ - - ]: 0 : if(Ma >= 0 && Ma < 1) { // Subsonic outflow
165 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
166 [ - - ]: 0 : auto gk = getDeformGrad(nmat, k, ul);
167 : 0 : ur[energyIdx(nmat, k)] = ul[volfracIdx(nmat, k)] *
168 : 0 : mat_blk[k].compute< EOS::totalenergy >(
169 [ - - ]: 0 : ur[densityIdx(nmat, k)]/ul[volfracIdx(nmat, k)], v1l, v2l, v3l, fp,
170 : : gk );
171 : : }
172 : :
173 : : // Internal cell primitive quantities using the separately reconstructed
174 : : // primitive quantities. This is used to get ghost state for primitive
175 : : // quantities
176 : :
177 : : // material pressures
178 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
179 : 0 : ur[ncomp+pressureIdx(nmat, k)] = ul[volfracIdx(nmat, k)] * fp;
180 : : }
181 : :
182 [ - - ][ - - ]: 0 : Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
[ - - ][ - - ]
183 : : "boundary state vector" );
184 : :
185 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
186 : : }
187 : :
188 : : //! \brief Boundary state function providing the left and right state of a
189 : : //! face at extrapolation boundaries
190 : : //! \param[in] ul Left (domain-internal) state
191 : : //! \return Left and right states for all scalar components in this PDE
192 : : //! system
193 : : //! \note The function signature must follow tk::StateFn. For multimat, the
194 : : //! left or right state is the vector of conserved quantities, followed by
195 : : //! the vector of primitive quantities appended to it.
196 : : static tk::StateFn::result_type
197 : 108016 : extrapolate( ncomp_t,
198 : : const std::vector< EOS >&,
199 : : const std::vector< tk::real >& ul,
200 : : tk::real, tk::real, tk::real, tk::real,
201 : : const std::array< tk::real, 3 >& )
202 : : {
203 : 108016 : return {{ ul, ul }};
204 : : }
205 : :
206 : : } // inciter::
207 : :
208 : : #endif // BCFunctions_h
|
