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