Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/MultiSpecies/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 "MiscMultiSpeciesFns.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.
29 : : static tk::StateFn::result_type
30 : 349200 : symmetry( [[maybe_unused]] ncomp_t ncomp,
31 : : const std::vector< EOS >&,
32 : : const std::vector< tk::real >& ul,
33 : : tk::real, tk::real, tk::real, tk::real,
34 : : const std::array< tk::real, 3 >& fn )
35 : : {
36 : 349200 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
37 : :
38 [ - + ][ - - ]: 349200 : Assert( ul.size() == ncomp, "Incorrect size for appended "
[ - - ][ - - ]
39 : : "internal state vector" );
40 : :
41 : 349200 : tk::real rho(0.0);
42 [ + + ]: 873000 : for (std::size_t k=0; k<nspec; ++k)
43 : 523800 : rho += ul[multispecies::densityIdx(nspec, k)];
44 : :
45 [ + - ]: 698400 : auto ur = ul;
46 : :
47 : : // Internal cell velocity components
48 : 349200 : auto v1l = ul[multispecies::momentumIdx(nspec, 0)]/rho;
49 : 349200 : auto v2l = ul[multispecies::momentumIdx(nspec, 1)]/rho;
50 : 349200 : auto v3l = ul[multispecies::momentumIdx(nspec, 2)]/rho;
51 : : // Normal component of velocity
52 : 349200 : auto vnl = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
53 : : // Ghost state velocity components
54 : 349200 : auto v1r = v1l - 2.0*vnl*fn[0];
55 : 349200 : auto v2r = v2l - 2.0*vnl*fn[1];
56 : 349200 : auto v3r = v3l - 2.0*vnl*fn[2];
57 : : // Boundary condition
58 : 349200 : ur[multispecies::momentumIdx(nspec, 0)] = rho * v1r;
59 : 349200 : ur[multispecies::momentumIdx(nspec, 1)] = rho * v2r;
60 : 349200 : ur[multispecies::momentumIdx(nspec, 2)] = rho * v3r;
61 : :
62 [ - + ][ - - ]: 349200 : Assert( ur.size() == ncomp, "Incorrect size for appended "
[ - - ][ - - ]
63 : : "boundary state vector" );
64 : :
65 [ + - ]: 698400 : return {{ std::move(ul), std::move(ur) }};
66 : : }
67 : :
68 : : //! \brief Boundary state function providing the left and right state of a
69 : : //! face at farfield boundaries
70 : : //! \param[in] ncomp Number of scalar components in this PDE system
71 : : //! \param[in] ul Left (domain-internal) state
72 : : //! \param[in] fn Unit face normal
73 : : //! \return Left and right states for all scalar components in this PDE
74 : : //! system
75 : : //! \details The farfield boudary calculation, implemented here, is
76 : : //! based on the characteristic theory of hyperbolic systems.
77 : : //! \note The function signature must follow tk::StateFn
78 : : static tk::StateFn::result_type
79 : 0 : farfield( [[maybe_unused]] ncomp_t ncomp,
80 : : const std::vector< EOS >& mat_blk,
81 : : const std::vector< tk::real >& ul,
82 : : tk::real, tk::real, tk::real, tk::real,
83 : : const std::array< tk::real, 3 >& fn )
84 : : {
85 : 0 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
86 : :
87 : : // Farfield primitive quantities
88 : : auto fp =
89 : 0 : g_inputdeck.get< tag::bc >()[0].get< tag::pressure >();
90 : : auto ft =
91 : 0 : g_inputdeck.get< tag::bc >()[0].get< tag::temperature >();
92 : : auto fu =
93 [ - - ]: 0 : g_inputdeck.get< tag::bc >()[0].get< tag::velocity >();
94 : : auto fspec =
95 [ - - ]: 0 : g_inputdeck.get< tag::bc >()[0].get< tag::mass_fractions >();
96 : :
97 [ - - ][ - - ]: 0 : Assert( ul.size() == ncomp, "Incorrect size for appended "
[ - - ][ - - ]
98 : : "internal state vector" );
99 : :
100 [ - - ]: 0 : auto ur = ul;
101 : :
102 : 0 : tk::real rhol(0.0);
103 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k)
104 : 0 : rhol += ul[multispecies::densityIdx(nspec, k)];
105 : :
106 : : // Internal cell velocity components
107 : 0 : auto v1l = ul[multispecies::momentumIdx(nspec, 0)]/rhol;
108 : 0 : auto v2l = ul[multispecies::momentumIdx(nspec, 1)]/rhol;
109 : 0 : auto v3l = ul[multispecies::momentumIdx(nspec, 2)]/rhol;
110 : :
111 : : // Normal velocity
112 : 0 : auto vn = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
113 : :
114 : : // Acoustic speed
115 [ - - ]: 0 : auto pl = mat_blk[0].compute< EOS::pressure >( rhol, v1l, v2l, v3l,
116 : 0 : ul[multispecies::energyIdx(nspec, 0)] );
117 [ - - ]: 0 : auto a = mat_blk[0].compute< EOS::soundspeed >( rhol, pl );
118 : :
119 : : // Mach number
120 : 0 : auto Ma = vn / a;
121 : :
122 [ - - ]: 0 : if (Ma <= -1) { // Supersonic inflow
123 : : // For supersonic inflow, all the characteristics are from outside.
124 : : // Therefore, we calculate the ghost cell state using the primitive
125 : : // variables from outside.
126 : 0 : tk::real rhor(0.0);
127 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k) {
128 [ - - ]: 0 : auto rhok = mat_blk[0].compute< EOS::density >(fp, ft);
129 : 0 : ur[multispecies::densityIdx(nspec,k)] = fspec[k] * rhok;
130 : 0 : rhor += ur[multispecies::densityIdx(nspec,k)];
131 : : }
132 : 0 : ur[multispecies::energyIdx(nspec,0)] =
133 [ - - ]: 0 : mat_blk[0].compute< EOS::totalenergy >(rhor, fu[0], fu[1], fu[2], fp);
134 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
135 : 0 : ur[multispecies::momentumIdx(nspec,i)] = rhor * fu[i];
136 : : }
137 : :
138 [ - - ][ - - ]: 0 : } else if (Ma > -1 && Ma < 0) { // Subsonic inflow
139 : : // For subsonic inflow, there is 1 outgoing characteristic and 4
140 : : // incoming characteristics. Therefore, we calculate the ghost cell state
141 : : // by taking pressure from the internal cell and other quantities from
142 : : // the outside.
143 : 0 : tk::real rhor(0.0);
144 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k) {
145 [ - - ]: 0 : auto rhok = mat_blk[0].compute< EOS::density >(pl, ft);
146 : 0 : ur[multispecies::densityIdx(nspec,k)] = fspec[k] * rhok;
147 : 0 : rhor += ur[multispecies::densityIdx(nspec,k)];
148 : : }
149 : 0 : ur[multispecies::energyIdx(nspec,0)] =
150 [ - - ]: 0 : mat_blk[0].compute< EOS::totalenergy >(rhor, fu[0], fu[1], fu[2], pl);
151 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
152 : 0 : ur[multispecies::momentumIdx(nspec,i)] = rhor * fu[i];
153 : 0 : }
154 : :
155 [ - - ][ - - ]: 0 : } else if (Ma >= 0 && Ma < 1) { // Subsonic outflow
156 : : // For subsonic outflow, there is 1 incoming characteristic and 4
157 : : // outgoing characteristics. Therefore, we calculate the ghost cell state
158 : : // by taking pressure from the outside and other quantities from the
159 : : // internal cell.
160 : 0 : ur[multispecies::energyIdx(nspec,0)] =
161 [ - - ]: 0 : mat_blk[0].compute< EOS::totalenergy >( rhol, v1l, v2l, v3l, fp );
162 : : }
163 : : // Otherwise, for supersonic outflow, all the characteristics are from
164 : : // internal cell. Therefore, we calculate the ghost cell state using the
165 : : // conservative variables from internal cell (which is what ur is
166 : : // initialized to).
167 : :
168 [ - - ][ - - ]: 0 : Assert( ur.size() == ncomp, "Incorrect size for appended "
[ - - ][ - - ]
169 : : "boundary state vector" );
170 : :
171 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
172 : : }
173 : :
174 : : //! \brief Boundary state function providing the left and right state of a
175 : : //! face at extrapolation boundaries
176 : : //! \param[in] ul Left (domain-internal) state
177 : : //! \return Left and right states for all scalar components in this PDE
178 : : //! system
179 : : //! \note The function signature must follow tk::StateFn.
180 : : static tk::StateFn::result_type
181 : 6750 : extrapolate( ncomp_t,
182 : : const std::vector< EOS >&,
183 : : const std::vector< tk::real >& ul,
184 : : tk::real, tk::real, tk::real, tk::real,
185 : : const std::array< tk::real, 3 >& )
186 : : {
187 : 6750 : return {{ ul, ul }};
188 : : }
189 : :
190 : : //! \brief Boundary state function providing the left and right state of a
191 : : //! face at no-slip wall boundaries
192 : : //! \param[in] ncomp Number of scalar components in this PDE system
193 : : //! \param[in] ul Left (domain-internal) state
194 : : // //! \param[in] fn Unit face normal
195 : : //! \return Left and right states for all scalar components in this PDE
196 : : //! system
197 : : //! \note The function signature must follow tk::StateFn.
198 : : static tk::StateFn::result_type
199 : 0 : noslipwall( [[maybe_unused]] ncomp_t ncomp,
200 : : const std::vector< EOS >&,
201 : : const std::vector< tk::real >& ul,
202 : : tk::real, tk::real, tk::real, tk::real,
203 : : const std::array< tk::real, 3 >& /*fn*/ )
204 : : {
205 : 0 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
206 : :
207 [ - - ][ - - ]: 0 : Assert( ul.size() == ncomp, "Incorrect size for appended "
[ - - ][ - - ]
208 : : "internal state vector" );
209 : :
210 : 0 : tk::real rho(0.0);
211 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k)
212 : 0 : rho += ul[multispecies::densityIdx(nspec, k)];
213 : :
214 [ - - ]: 0 : auto ur = ul;
215 : :
216 : : // Internal cell velocity components
217 : 0 : auto v1l = ul[multispecies::momentumIdx(nspec, 0)]/rho;
218 : 0 : auto v2l = ul[multispecies::momentumIdx(nspec, 1)]/rho;
219 : 0 : auto v3l = ul[multispecies::momentumIdx(nspec, 2)]/rho;
220 : : // Ghost state velocity components
221 : 0 : auto v1r = -v1l;
222 : 0 : auto v2r = -v2l;
223 : 0 : auto v3r = -v3l;
224 : : // Boundary condition
225 : 0 : ur[multispecies::momentumIdx(nspec, 0)] = rho * v1r;
226 : 0 : ur[multispecies::momentumIdx(nspec, 1)] = rho * v2r;
227 : 0 : ur[multispecies::momentumIdx(nspec, 2)] = rho * v3r;
228 : :
229 [ - - ][ - - ]: 0 : Assert( ur.size() == ncomp, "Incorrect size for appended "
[ - - ][ - - ]
230 : : "boundary state vector" );
231 : :
232 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
233 : : }
234 : :
235 : : //----------------------------------------------------------------------------
236 : : // Boundary Gradient functions
237 : : //----------------------------------------------------------------------------
238 : :
239 : : //! \brief Boundary gradient function copying the left gradient to the right
240 : : //! gradient at a face
241 : : //! \param[in] dul Left (domain-internal) state
242 : : //! \return Left and right states for all scalar components in this PDE
243 : : //! system
244 : : //! \note The function signature must follow tk::StateFn.
245 : : static tk::StateFn::result_type
246 : 0 : noOpGrad( ncomp_t,
247 : : const std::vector< EOS >&,
248 : : const std::vector< tk::real >& dul,
249 : : tk::real, tk::real, tk::real, tk::real,
250 : : const std::array< tk::real, 3 >& )
251 : : {
252 : 0 : return {{ dul, dul }};
253 : : }
254 : :
255 : : //! \brief Boundary gradient function for the symmetry boundary condition
256 : : //! \param[in] ncomp Number of variables whos gradients are needed
257 : : //! \param[in] dul Left (domain-internal) gradients
258 : : //! \return Left and right states for all scalar components in this PDE
259 : : //! system
260 : : //! \note The function signature must follow tk::StateFn.
261 : : static tk::StateFn::result_type
262 : 0 : symmetryGrad( ncomp_t ncomp,
263 : : const std::vector< EOS >&,
264 : : const std::vector< tk::real >& dul,
265 : : tk::real, tk::real, tk::real, tk::real,
266 : : const std::array< tk::real, 3 >& )
267 : : {
268 [ - - ][ - - ]: 0 : Assert(dul.size() == 3*ncomp, "Incorrect size of boundary gradient vector");
[ - - ][ - - ]
269 : :
270 [ - - ]: 0 : auto dur = dul;
271 : :
272 [ - - ]: 0 : for (std::size_t i=0; i<3*ncomp; ++i)
273 : 0 : dur[i] = -dul[i];
274 : :
275 [ - - ]: 0 : return {{ std::move(dul), std::move(dur) }};
276 : : }
277 : :
278 : : } // inciter::
279 : :
280 : : #endif // BCFunctions_h
|
