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 : : Assert( ul.size() == ncomp, "Incorrect size for appended "
39 : : "internal state vector" );
40 : :
41 : : 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 : 349200 : 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 : : Assert( ur.size() == ncomp, "Incorrect size for appended "
63 : : "boundary state vector" );
64 : :
65 [ + - ]: 349200 : 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 : : Assert( ul.size() == ncomp, "Incorrect size for appended "
98 : : "internal state vector" );
99 : :
100 [ - - ]: 0 : auto ur = ul;
101 : :
102 [ - - ]: 0 : Mixture mixl(nspec, ul, mat_blk); // Initialize mixture class
103 : 0 : tk::real rhol = mixl.get_mix_density();
104 : :
105 : : // Internal cell velocity components
106 [ - - ]: 0 : auto v1l = ul[multispecies::momentumIdx(nspec, 0)]/rhol;
107 : 0 : auto v2l = ul[multispecies::momentumIdx(nspec, 1)]/rhol;
108 : 0 : auto v3l = ul[multispecies::momentumIdx(nspec, 2)]/rhol;
109 : :
110 : : // Normal velocity
111 [ - - ]: 0 : auto vn = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
112 : :
113 : : // Acoustic speed
114 [ - - ]: 0 : auto pl = mixl.pressure(rhol, v1l, v2l, v3l,
115 : : ul[multispecies::energyIdx(nspec,0)], mat_blk);
116 [ - - ]: 0 : auto a = mixl.frozen_soundspeed(rhol, pl, mat_blk);
117 : :
118 : : // Mach number
119 : 0 : auto Ma = vn / a;
120 : :
121 [ - - ]: 0 : if (Ma <= -1) { // Supersonic inflow
122 : : // For supersonic inflow, all the characteristics are from outside.
123 : : // Therefore, we calculate the ghost cell state using the primitive
124 : : // variables from outside.
125 [ - - ]: 0 : Mixture mixr(nspec, fspec, fp, ft, mat_blk);
126 : :
127 : 0 : tk::real rhor = mixr.get_mix_density();
128 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k) {
129 : 0 : ur[multispecies::densityIdx(nspec,k)] = fspec[k] * rhor;
130 : : }
131 : 0 : ur[multispecies::energyIdx(nspec,0)] = mixr.totalenergy(rhor, fu[0],
132 [ - - ]: 0 : fu[1], fu[2], fp, mat_blk);
133 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
134 : 0 : ur[multispecies::momentumIdx(nspec,i)] = rhor * fu[i];
135 : : }
136 : :
137 [ - - ][ - - ]: 0 : } else if (Ma > -1 && Ma < 0) { // Subsonic inflow
138 : : // For subsonic inflow, there is 1 outgoing characteristic and 4
139 : : // incoming characteristics. Therefore, we calculate the ghost cell state
140 : : // by taking pressure from the internal cell and other quantities from
141 : : // the outside.
142 [ - - ]: 0 : Mixture mixr(nspec, fspec, pl, ft, mat_blk);
143 : :
144 : 0 : tk::real rhor = mixr.get_mix_density();
145 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k) {
146 : 0 : ur[multispecies::densityIdx(nspec,k)] = fspec[k] * rhor;
147 : : }
148 : 0 : ur[multispecies::energyIdx(nspec,0)] = mixr.totalenergy(rhor, fu[0],
149 [ - - ]: 0 : fu[1], fu[2], pl, mat_blk);
150 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
151 : 0 : ur[multispecies::momentumIdx(nspec,i)] = rhor * fu[i];
152 : 0 : }
153 : :
154 [ - - ][ - - ]: 0 : } else if (Ma >= 0 && Ma < 1) { // Subsonic outflow
155 : : // For subsonic outflow, there is 1 incoming characteristic and 4
156 : : // outgoing characteristics. Therefore, we calculate the ghost cell state
157 : : // by taking pressure from the outside and other quantities from the
158 : : // internal cell.
159 : : std::vector< tk::real > massfrac_l;
160 [ - - ]: 0 : for (std::size_t k = 0; k < nspec; k++)
161 : 0 : massfrac_l[k] = ul[multispecies::densityIdx(nspec, k)] / rhol;
162 [ - - ]: 0 : tk::real tl = mixl.temperature(rhol, v1l, v2l, v3l,
163 [ - - ]: 0 : ul[multispecies::energyIdx(nspec,0)], mat_blk);
164 [ - - ]: 0 : Mixture mixr(nspec, massfrac_l, fp, tl, mat_blk);
165 : :
166 [ - - ][ - - ]: 0 : ur[multispecies::energyIdx(nspec,0)] = mixr.totalenergy(rhol, v1l,
167 : : v2l, v3l, fp, mat_blk);
168 : : }
169 : : // Otherwise, for supersonic outflow, all the characteristics are from
170 : : // internal cell. Therefore, we calculate the ghost cell state using the
171 : : // conservative variables from internal cell (which is what ur is
172 : : // initialized to).
173 : :
174 : : Assert( ur.size() == ncomp, "Incorrect size for appended "
175 : : "boundary state vector" );
176 : :
177 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
178 : : }
179 : :
180 : : //! \brief Boundary state function providing the left and right state of a
181 : : //! face at extrapolation boundaries
182 : : //! \param[in] ul Left (domain-internal) state
183 : : //! \return Left and right states for all scalar components in this PDE
184 : : //! system
185 : : //! \note The function signature must follow tk::StateFn.
186 : : static tk::StateFn::result_type
187 : 6750 : extrapolate( ncomp_t,
188 : : const std::vector< EOS >&,
189 : : const std::vector< tk::real >& ul,
190 : : tk::real, tk::real, tk::real, tk::real,
191 : : const std::array< tk::real, 3 >& )
192 : : {
193 : 6750 : return {{ ul, ul }};
194 : : }
195 : :
196 : : //! \brief Boundary state function providing the left and right state of a
197 : : //! face at no-slip wall boundaries
198 : : //! \param[in] ncomp Number of scalar components in this PDE system
199 : : //! \param[in] ul Left (domain-internal) state
200 : : // //! \param[in] fn Unit face normal
201 : : //! \return Left and right states for all scalar components in this PDE
202 : : //! system
203 : : //! \note The function signature must follow tk::StateFn.
204 : : static tk::StateFn::result_type
205 : 0 : noslipwall( [[maybe_unused]] ncomp_t ncomp,
206 : : const std::vector< EOS >&,
207 : : const std::vector< tk::real >& ul,
208 : : tk::real, tk::real, tk::real, tk::real,
209 : : const std::array< tk::real, 3 >& /*fn*/ )
210 : : {
211 : 0 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
212 : :
213 : : Assert( ul.size() == ncomp, "Incorrect size for appended "
214 : : "internal state vector" );
215 : :
216 : : tk::real rho(0.0);
217 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k)
218 : 0 : rho += ul[multispecies::densityIdx(nspec, k)];
219 : :
220 : 0 : auto ur = ul;
221 : :
222 : : // Internal cell velocity components
223 [ - - ]: 0 : auto v1l = ul[multispecies::momentumIdx(nspec, 0)]/rho;
224 : 0 : auto v2l = ul[multispecies::momentumIdx(nspec, 1)]/rho;
225 : 0 : auto v3l = ul[multispecies::momentumIdx(nspec, 2)]/rho;
226 : : // Ghost state velocity components
227 : 0 : auto v1r = -v1l;
228 : 0 : auto v2r = -v2l;
229 : 0 : auto v3r = -v3l;
230 : : // Boundary condition
231 [ - - ]: 0 : ur[multispecies::momentumIdx(nspec, 0)] = rho * v1r;
232 : 0 : ur[multispecies::momentumIdx(nspec, 1)] = rho * v2r;
233 : 0 : ur[multispecies::momentumIdx(nspec, 2)] = rho * v3r;
234 : :
235 : : Assert( ur.size() == ncomp, "Incorrect size for appended "
236 : : "boundary state vector" );
237 : :
238 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
239 : : }
240 : :
241 : : //----------------------------------------------------------------------------
242 : : // Boundary Gradient functions
243 : : //----------------------------------------------------------------------------
244 : :
245 : : //! \brief Boundary gradient function copying the left gradient to the right
246 : : //! gradient at a face
247 : : //! \param[in] dul Left (domain-internal) state
248 : : //! \return Left and right states for all scalar components in this PDE
249 : : //! system
250 : : //! \note The function signature must follow tk::StateFn.
251 : : static tk::StateFn::result_type
252 : 0 : noOpGrad( ncomp_t,
253 : : const std::vector< EOS >&,
254 : : const std::vector< tk::real >& dul,
255 : : tk::real, tk::real, tk::real, tk::real,
256 : : const std::array< tk::real, 3 >& )
257 : : {
258 : 0 : return {{ dul, dul }};
259 : : }
260 : :
261 : : //! \brief Boundary gradient function for the symmetry boundary condition
262 : : //! \param[in] ncomp Number of variables whos gradients are needed
263 : : //! \param[in] dul Left (domain-internal) gradients
264 : : //! \return Left and right states for all scalar components in this PDE
265 : : //! system
266 : : //! \note The function signature must follow tk::StateFn.
267 : : static tk::StateFn::result_type
268 : 0 : symmetryGrad( ncomp_t ncomp,
269 : : const std::vector< EOS >&,
270 : : const std::vector< tk::real >& dul,
271 : : tk::real, tk::real, tk::real, tk::real,
272 : : const std::array< tk::real, 3 >& )
273 : : {
274 : : Assert(dul.size() == 3*ncomp, "Incorrect size of boundary gradient vector");
275 : :
276 : 0 : auto dur = dul;
277 : :
278 [ - - ]: 0 : for (std::size_t i=0; i<3*ncomp; ++i)
279 : 0 : dur[i] = -dul[i];
280 : :
281 [ - - ]: 0 : return {{ std::move(dul), std::move(dur) }};
282 : : }
283 : :
284 : : } // inciter::
285 : :
286 : : #endif // BCFunctions_h
|
