Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Physics/BC.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.,
7 : : 2022-2024 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief Boundary conditions
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "BC.hpp"
14 : : #include "EOS.hpp"
15 : : #include "Box.hpp"
16 : : #include "Problems.hpp"
17 : : #include "InciterConfig.hpp"
18 : :
19 : : namespace inciter {
20 : :
21 : : extern ctr::Config g_cfg;
22 : :
23 : : } // ::inciter
24 : :
25 : : namespace physics {
26 : :
27 : : using inciter::g_cfg;
28 : :
29 : : void
30 [ + + ]: 77800 : dirbc( tk::Fields& U,
31 : : tk::real t,
32 : : const std::array< std::vector< tk::real >, 3 >& coord,
33 : : const std::vector< std::unordered_set< std::size_t > >& boxnodes,
34 : : const std::vector< std::size_t >& dirbcmask,
35 : : const std::vector< double >& dirbcval )
36 : : // *****************************************************************************
37 : : // Set Dirichlet boundary conditions at nodes
38 : : //! \param[in] t Physical time at which to evaluate BCs
39 : : //! \param[in] U Solution vector at recent time step
40 : : //! \param[in] coord Mesh node coordinates
41 : : //! \param[in] boxnodes List of nodes at which box user ICs are set
42 : : //! \param[in] dirbcmask Nodes and component masks for Dirichlet BCs
43 : : //! \param[in] dirbcval Nodes and component values for Dirichlet BCs
44 : : // *****************************************************************************
45 : : {
46 [ + + ]: 77800 : if (g_cfg.get< tag::bc_dir >().empty()) return;
47 : :
48 : : auto ncomp = U.nprop();
49 : 34650 : auto nmask = ncomp + 1;
50 : :
51 : : Assert( dirbcmask.size() % nmask == 0, "Size mismatch" );
52 : : Assert( dirbcval.size() % nmask == 0, "Size mismatch" );
53 : :
54 : : const auto& x = coord[0];
55 : : const auto& y = coord[1];
56 : : const auto& z = coord[2];
57 : 34650 : auto ic = problems::IC();
58 : :
59 [ + + ]: 2897896 : for (std::size_t i=0; i<dirbcmask.size()/nmask; ++i) {
60 [ - + ]: 2863246 : auto p = dirbcmask[i*nmask+0]; // local node id
61 [ - + ]: 2863246 : auto u = ic( x[p], y[p], z[p], t ); // evaluate solution/ic
62 [ + - ]: 2863246 : problems::box( p, u, boxnodes ); // overwrite with box value
63 [ + + ]: 16095294 : for (std::size_t c=0; c<ncomp; ++c) {
64 [ + + ]: 13232048 : auto mask = dirbcmask[i*nmask+1+c];
65 [ + + ]: 13232048 : if (mask == 1) { // mask == 1: IC+box value
66 : 10970411 : U(p,c) = u[c];
67 [ + + ][ + - ]: 2261637 : } else if (mask == 2 && !dirbcval.empty()) { // mask == 2: BC value
68 : 27114 : U(p,c) = dirbcval[i*nmask+1+c];
69 : : }
70 : : }
71 : : }
72 : : }
73 : :
74 : : void
75 [ + + ]: 4660 : dirbcp( tk::Fields& U,
76 : : const std::array< std::vector< tk::real >, 3 >& coord,
77 : : const std::vector< std::size_t >& dirbcmaskp,
78 : : const std::vector< double >& dirbcvalp )
79 : : // *****************************************************************************
80 : : // Set pressure Dirichlet boundary conditions at nodes
81 : : //! \param[in] U Solution vector at recent time step
82 : : //! \param[in] coord Mesh node coordinates
83 : : //! \param[in] dirbcmaskp Nodes and component masks for Dirichlet BCs
84 : : //! \param[in] dirbcvalp Nodes and component values for Dirichlet BCs
85 : : // *****************************************************************************
86 : : {
87 [ + + ]: 4660 : if (g_cfg.get< tag::pre_bc_dir >().empty()) return;
88 : :
89 : : std::size_t nmask = 1 + 1;
90 : :
91 : : Assert( dirbcmaskp.size() % nmask == 0, "Size mismatch" );
92 : : Assert( dirbcvalp.size() % nmask == 0, "Size mismatch" );
93 : :
94 : : const auto& x = coord[0];
95 : : const auto& y = coord[1];
96 : : const auto& z = coord[2];
97 : 1520 : auto ic = problems::PRESSURE_IC();
98 : :
99 [ + + ]: 10640 : for (std::size_t i=0; i<dirbcmaskp.size()/nmask; ++i) {
100 [ - + ]: 9120 : auto p = dirbcmaskp[i*nmask+0]; // local node id
101 : 9120 : auto mask = dirbcmaskp[i*nmask+1];
102 [ - + ]: 9120 : if (mask == 1) { // mask == 1: IC value
103 [ - - ]: 0 : U(p,0) = ic( x[p], y[p], z[p] );
104 [ + - ][ + - ]: 9120 : } else if (mask == 2 && !dirbcvalp.empty()) { // mask == 2: BC value
105 : 9120 : U(p,0) = dirbcvalp[i*nmask+1];
106 : : }
107 : : }
108 : : }
109 : :
110 : : void
111 [ + + ]: 78458 : symbc( tk::Fields& U,
112 : : const std::vector< std::size_t >& symbcnodes,
113 : : const std::vector< tk::real >& symbcnorms,
114 : : std::size_t pos )
115 : : // *****************************************************************************
116 : : // Set symmetry boundary conditions at nodes
117 : : //! \param[in] U Solution vector at recent time step
118 : : //! \param[in] symbcnodes Node ids at which to set symmetry BCs
119 : : //! \param[in] symbcnorms Normals at nodes at which to set symmetry BCs
120 : : //! \param[in] pos Position at which the three velocity components are in U
121 : : // *****************************************************************************
122 : : {
123 [ + + ]: 78458 : if (g_cfg.get< tag::bc_sym >().empty()) return;
124 : :
125 : : Assert( symbcnodes.size()*3 == symbcnorms.size(), "Size mismatch" );
126 : :
127 [ + + ]: 2383980 : for (std::size_t i=0; i<symbcnodes.size(); ++i) {
128 : 2339623 : auto p = symbcnodes[i];
129 : 2339623 : auto n = symbcnorms.data() + i*3;
130 : : auto& u = U(p,pos+0);
131 : 2339623 : auto& v = U(p,pos+1);
132 : 2339623 : auto& w = U(p,pos+2);
133 : 2339623 : auto vn = u*n[0] + v*n[1] + w*n[2];
134 : 2339623 : u -= vn * n[0];
135 : 2339623 : v -= vn * n[1];
136 : 2339623 : w -= vn * n[2];
137 : : }
138 : : }
139 : :
140 : : void
141 [ + + ]: 13420 : noslipbc( tk::Fields& U,
142 : : const std::vector< std::size_t >& noslipbcnodes,
143 : : std::size_t pos )
144 : : // *****************************************************************************
145 : : // Set noslip boundary conditions at nodes
146 : : //! \param[in] U Solution vector at recent time step
147 : : //! \param[in] noslipbcnodes Node ids at which to set noslip BCs
148 : : //! \param[in] pos Position at which the three velocity components are in U
149 : : // *****************************************************************************
150 : : {
151 [ + + ]: 13420 : if (g_cfg.get< tag::bc_noslip >().empty()) return;
152 : :
153 [ + + ]: 426322 : for (auto p : noslipbcnodes) U(p,pos+0) = U(p,pos+1) = U(p,pos+2) = 0.0;
154 : : }
155 : :
156 : : void
157 [ + + ]: 64380 : farbc( tk::Fields& U,
158 : : const std::vector< std::size_t >& farbcnodes,
159 : : const std::vector< tk::real >& farbcnorms )
160 : : // *****************************************************************************
161 : : // Set farfield boundary conditions at nodes
162 : : //! \param[in] U Solution vector at recent time step
163 : : //! \param[in] farbcnodes Nodes ids at which to set farfield BCs
164 : : //! \param[in] farbcnorms Normals at nodes at which to set farfield BCs
165 : : // *****************************************************************************
166 : : {
167 [ + + ]: 64380 : if (g_cfg.get< tag::bc_far >().empty()) return;
168 : :
169 : : Assert( farbcnodes.size()*3 == farbcnorms.size(), "Size mismatch" );
170 : :
171 : : // cppcheck-suppress unreadVariable
172 [ - + ]: 2270 : tk::real fr = g_cfg.get< tag::bc_far_density >();
173 : :
174 : : const auto& fue = g_cfg.get< tag::bc_far_velocity >();
175 [ - + ][ - - ]: 2270 : ErrChk( !fue.empty(), "No farfield velocity specified" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
176 : : // cppcheck-suppress unreadVariable
177 : 2270 : tk::real fu = fue[0];
178 : : // cppcheck-suppress unreadVariable
179 : 2270 : tk::real fv = fue[1];
180 : : // cppcheck-suppress unreadVariable
181 : 2270 : tk::real fw = fue[2];
182 : :
183 : 2270 : tk::real fp = g_cfg.get< tag::bc_far_pressure >();
184 : :
185 [ + + ]: 45889 : for (std::size_t i=0; i<farbcnodes.size(); ++i) {
186 : 43619 : auto p = farbcnodes[i];
187 : 43619 : auto nx = farbcnorms[i*3+0];
188 : 43619 : auto ny = farbcnorms[i*3+1];
189 [ - + ]: 43619 : auto nz = farbcnorms[i*3+2];
190 : : auto& r = U(p,0);
191 : : auto& ru = U(p,1);
192 : : auto& rv = U(p,2);
193 : : auto& rw = U(p,3);
194 : : auto& re = U(p,4);
195 : : //auto vn = (ru*nx + rv*ny + rw*nz)/r;
196 [ - + ]: 43619 : auto vn = fu*nx + fv*ny + fw*nz;
197 : : //auto a = eos::soundspeed( r,
198 : : // eos::pressure( re - 0.5*(ru*ru + rv*rv + rw*rw)/r ) );
199 : : auto a = eos::soundspeed( fr, fp );
200 : 43619 : auto M = vn / a;
201 [ - + ]: 43619 : if (M <= -1.0) {
202 : : // supersonic inflow, all characteristics from outside
203 : 0 : r = fr;
204 : 0 : ru = fr * fu;
205 : 0 : rv = fr * fv;
206 : 0 : rw = fr * fw;
207 : 0 : re = eos::totalenergy( fr, fu, fv, fw, fp );
208 [ + - ][ + + ]: 43619 : } else if (M > -1.0 && M < 0.0) {
209 : : // subsonic inflow: 1 outgoing and 4 incoming characteristics,
210 : : // pressure from inside, rest from outside
211 : 21799 : auto pr = eos::pressure( re - 0.5*(ru*ru + rv*rv + rw*rw)/r );
212 : 21799 : r = fr;
213 : 21799 : ru = fr * fu;
214 : 21799 : rv = fr * fv;
215 : 21799 : rw = fr * fw;
216 : 21799 : re = eos::totalenergy( fr, fu, fv, fw, pr );
217 [ + - ][ + - ]: 43619 : } else if (M >= 0.0 && M < 1.0) {
218 : : // subsonic outflow: 1 incoming and 4 outgoing characteristics,
219 : : // pressure from outside, rest from inside
220 : 21820 : re = eos::totalenergy( r, ru/r, rv/r, rw/r, fp );
221 : : }
222 : : }
223 : : }
224 : :
225 : : void
226 : 64380 : prebc( tk::Fields& U,
227 : : const std::vector< std::size_t >& prebcnodes,
228 : : const std::vector< tk::real >& prebcvals )
229 : : // *****************************************************************************
230 : : // Set pressure boundary conditions at nodes
231 : : //! \param[in] U Solution vector at recent time step
232 : : //! \param[in] prebcnodes Node ids at which to set pressure BCs
233 : : //! \param[in] prebcvals Density and pressure values at pressure BC nodes
234 : : // *****************************************************************************
235 : : {
236 : : Assert( prebcnodes.size()*2 == prebcvals.size(), "Size mismatch" );
237 : :
238 [ + + ]: 74714 : for (std::size_t i=0; i<prebcnodes.size(); ++i) {
239 : 10334 : auto p = prebcnodes[i];
240 : 10334 : U(p,0) = prebcvals[i*2+0];
241 : 10334 : U(p,4) = eos::totalenergy( U(p,0), U(p,1)/U(p,0), U(p,2)/U(p,0),
242 : 10334 : U(p,3)/U(p,0), prebcvals[i*2+1] );
243 : : }
244 : 64380 : }
245 : :
246 : : } // physics::
|
