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+1, "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+1, "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+1, "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 Tl = ul[ncomp+multispecies::temperatureIdx(nspec,0)];
115 [ - - ]: 0 : auto pl = mixl.pressure(rhol, Tl);
116 [ - - ]: 0 : auto a = mixl.frozen_soundspeed(rhol, Tl, 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], ft, 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], Tl, 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 : Mixture mixr(nspec, massfrac_l, fp, Tl, mat_blk);
163 : :
164 [ - - ][ - - ]: 0 : ur[multispecies::energyIdx(nspec,0)] = mixr.totalenergy(rhol, v1l,
165 : : v2l, v3l, fp, mat_blk);
166 : : }
167 : : // Otherwise, for supersonic outflow, all the characteristics are from
168 : : // internal cell. Therefore, we calculate the ghost cell state using the
169 : : // conservative variables from internal cell (which is what ur is
170 : : // initialized to).
171 : :
172 : : Assert( ur.size() == ncomp+1, "Incorrect size for appended "
173 : : "boundary state vector" );
174 : :
175 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
176 : : }
177 : :
178 : : //! \brief Boundary state function providing the left and right state of a
179 : : //! face at extrapolation boundaries
180 : : //! \param[in] ul Left (domain-internal) state
181 : : //! \return Left and right states for all scalar components in this PDE
182 : : //! system
183 : : //! \note The function signature must follow tk::StateFn.
184 : : static tk::StateFn::result_type
185 : 6750 : extrapolate( ncomp_t,
186 : : const std::vector< EOS >&,
187 : : const std::vector< tk::real >& ul,
188 : : tk::real, tk::real, tk::real, tk::real,
189 : : const std::array< tk::real, 3 >& )
190 : : {
191 : 6750 : return {{ ul, ul }};
192 : : }
193 : :
194 : : //! \brief Boundary state function providing the left and right state of a
195 : : //! face at no-slip wall boundaries
196 : : //! \param[in] ncomp Number of scalar components in this PDE system
197 : : //! \param[in] ul Left (domain-internal) state
198 : : // //! \param[in] fn Unit face normal
199 : : //! \return Left and right states for all scalar components in this PDE
200 : : //! system
201 : : //! \note The function signature must follow tk::StateFn.
202 : : static tk::StateFn::result_type
203 : 0 : noslipwall( [[maybe_unused]] ncomp_t ncomp,
204 : : const std::vector< EOS >&,
205 : : const std::vector< tk::real >& ul,
206 : : tk::real, tk::real, tk::real, tk::real,
207 : : const std::array< tk::real, 3 >& /*fn*/ )
208 : : {
209 : 0 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
210 : :
211 : : Assert( ul.size() == ncomp+1, "Incorrect size for appended "
212 : : "internal state vector" );
213 : :
214 : : tk::real rho(0.0);
215 [ - - ]: 0 : for (std::size_t k=0; k<nspec; ++k)
216 : 0 : rho += ul[multispecies::densityIdx(nspec, k)];
217 : :
218 : 0 : auto ur = ul;
219 : :
220 : : // Internal cell velocity components
221 [ - - ]: 0 : auto v1l = ul[multispecies::momentumIdx(nspec, 0)]/rho;
222 : 0 : auto v2l = ul[multispecies::momentumIdx(nspec, 1)]/rho;
223 : 0 : auto v3l = ul[multispecies::momentumIdx(nspec, 2)]/rho;
224 : : // Ghost state velocity components
225 : 0 : auto v1r = -v1l;
226 : 0 : auto v2r = -v2l;
227 : 0 : auto v3r = -v3l;
228 : : // Boundary condition
229 [ - - ]: 0 : ur[multispecies::momentumIdx(nspec, 0)] = rho * v1r;
230 : 0 : ur[multispecies::momentumIdx(nspec, 1)] = rho * v2r;
231 : 0 : ur[multispecies::momentumIdx(nspec, 2)] = rho * v3r;
232 : :
233 : : Assert( ur.size() == ncomp+1, "Incorrect size for appended "
234 : : "boundary state vector" );
235 : :
236 [ - - ]: 0 : return {{ std::move(ul), std::move(ur) }};
237 : : }
238 : :
239 : : //----------------------------------------------------------------------------
240 : : // Boundary Gradient functions
241 : : //----------------------------------------------------------------------------
242 : :
243 : : //! \brief Boundary gradient function copying the left gradient to the right
244 : : //! gradient at a face
245 : : //! \param[in] dul Left (domain-internal) state
246 : : //! \return Left and right states for all scalar components in this PDE
247 : : //! system
248 : : //! \note The function signature must follow tk::StateFn.
249 : : static tk::StateFn::result_type
250 : 0 : noOpGrad( ncomp_t,
251 : : const std::vector< EOS >&,
252 : : const std::vector< tk::real >& dul,
253 : : tk::real, tk::real, tk::real, tk::real,
254 : : const std::array< tk::real, 3 >& )
255 : : {
256 : 0 : return {{ dul, dul }};
257 : : }
258 : :
259 : : //! \brief Boundary gradient function for the symmetry boundary condition
260 : : //! \param[in] ncomp Number of variables whos gradients are needed
261 : : //! \param[in] dul Left (domain-internal) gradients
262 : : //! \return Left and right states for all scalar components in this PDE
263 : : //! system
264 : : //! \note The function signature must follow tk::StateFn.
265 : : static tk::StateFn::result_type
266 : 0 : symmetryGrad( ncomp_t ncomp,
267 : : const std::vector< EOS >&,
268 : : const std::vector< tk::real >& dul,
269 : : tk::real, tk::real, tk::real, tk::real,
270 : : const std::array< tk::real, 3 >& )
271 : : {
272 : : Assert(dul.size() == 3*ncomp, "Incorrect size of boundary gradient vector");
273 : :
274 : 0 : auto dur = dul;
275 : :
276 [ - - ]: 0 : for (std::size_t i=0; i<3*ncomp; ++i)
277 : 0 : dur[i] = -dul[i];
278 : :
279 [ - - ]: 0 : return {{ std::move(dul), std::move(dur) }};
280 : : }
281 : :
282 : : } // inciter::
283 : :
284 : : #endif // BCFunctions_h
|
