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-2025 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 : 94904 : dirbc( std::size_t meshid,
31 : : tk::Fields& U,
32 : : tk::real t,
33 : : const std::array< std::vector< tk::real >, 3 >& coord,
34 : : const std::vector< std::unordered_set< std::size_t > >& boxnodes,
35 : : const std::vector< std::size_t >& dirbcmask,
36 : : const std::vector< double >& dirbcval )
37 : : // *****************************************************************************
38 : : // Set Dirichlet boundary conditions at nodes
39 : : //! \param[in] meshid Mesh id to use
40 : : //! \param[in] U Solution vector at recent time step
41 : : //! \param[in] t Physical time at which to evaluate BCs
42 : : //! \param[in] coord Mesh node coordinates
43 : : //! \param[in] boxnodes List of nodes at which box user ICs are set
44 : : //! \param[in] dirbcmask Nodes and component masks for Dirichlet BCs
45 : : //! \param[in] dirbcval Nodes and component values for Dirichlet BCs
46 : : // *****************************************************************************
47 : : {
48 : : auto ncomp = U.nprop();
49 : 94904 : 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 : 94904 : auto ic = problems::IC();
58 : :
59 [ + + ]: 4702402 : for (std::size_t i=0; i<dirbcmask.size()/nmask; ++i) {
60 [ - + ]: 4607498 : auto p = dirbcmask[i*nmask+0]; // local node id
61 [ - + ]: 4607498 : auto u = ic( x[p], y[p], z[p], t, meshid ); // evaluate solution/ic
62 [ + - ]: 4607498 : problems::box( p, u, boxnodes ); // overwrite with box value
63 [ + + ]: 25599466 : for (std::size_t c=0; c<ncomp; ++c) {
64 [ + + ]: 20991968 : auto mask = dirbcmask[i*nmask+1+c];
65 [ + + ]: 20991968 : if (mask == 1) { // mask == 1: IC+box value
66 : 16037967 : U(p,c) = u[c];
67 [ + + ][ + - ]: 4954001 : } else if (mask == 2 && !dirbcval.empty()) { // mask == 2: BC value
68 : 55374 : U(p,c) = dirbcval[i*nmask+1+c];
69 : : }
70 : : }
71 : : }
72 : 94904 : }
73 : :
74 : : void
75 : 17200 : dirbcp( std::size_t meshid,
76 : : tk::Fields& U,
77 : : const std::array< std::vector< tk::real >, 3 >& coord,
78 : : const std::vector< std::size_t >& dirbcmaskp,
79 : : const std::vector< double >& dirbcvalp )
80 : : // *****************************************************************************
81 : : // Set pressure Dirichlet boundary conditions at nodes
82 : : //! \param[in] meshid Mesh id to use
83 : : //! \param[in] U Solution vector at recent time step
84 : : //! \param[in] coord Mesh node coordinates
85 : : //! \param[in] dirbcmaskp Nodes and component masks for Dirichlet BCs
86 : : //! \param[in] dirbcvalp Nodes and component values for Dirichlet BCs
87 : : // *****************************************************************************
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 : 17200 : auto ic = problems::PRESSURE_IC();
98 : :
99 [ + + ]: 47920 : for (std::size_t i=0; i<dirbcmaskp.size()/nmask; ++i) {
100 : 30720 : auto p = dirbcmaskp[i*nmask+0]; // local node id
101 : 30720 : auto mask = dirbcmaskp[i*nmask+1];
102 [ + + ]: 30720 : if (mask == 1) { // mask == 1: IC value
103 [ - + ]: 30720 : U(p,0) = ic( x[p], y[p], z[p], meshid );
104 [ + - ][ + - ]: 15360 : } else if (mask == 2 && !dirbcvalp.empty()) { // mask == 2: BC value
105 : 15360 : U(p,0) = dirbcvalp[i*nmask+1];
106 : : }
107 : : }
108 : 17200 : }
109 : :
110 : : void
111 : 95674 : 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 : : Assert( symbcnodes.size()*3 == symbcnorms.size(), "Size mismatch" );
124 : :
125 [ + + ]: 2664715 : for (std::size_t i=0; i<symbcnodes.size(); ++i) {
126 : 2569041 : auto p = symbcnodes[i];
127 : 2569041 : auto n = symbcnorms.data() + i*3;
128 : : auto& u = U(p,pos+0);
129 : 2569041 : auto& v = U(p,pos+1);
130 : 2569041 : auto& w = U(p,pos+2);
131 : 2569041 : auto vn = u*n[0] + v*n[1] + w*n[2];
132 : 2569041 : u -= vn * n[0];
133 : 2569041 : v -= vn * n[1];
134 : 2569041 : w -= vn * n[2];
135 : : }
136 : 95674 : }
137 : :
138 : : void
139 : 30404 : noslipbc( tk::Fields& U,
140 : : const std::vector< std::size_t >& noslipbcnodes,
141 : : std::size_t pos )
142 : : // *****************************************************************************
143 : : // Set noslip boundary conditions at nodes
144 : : //! \param[in] U Solution vector at recent time step
145 : : //! \param[in] noslipbcnodes Node ids at which to set noslip BCs
146 : : //! \param[in] pos Position at which the three velocity components are in U
147 : : // *****************************************************************************
148 : : {
149 [ + + ]: 783674 : for (auto p : noslipbcnodes) U(p,pos+0) = U(p,pos+1) = U(p,pos+2) = 0.0;
150 : 30404 : }
151 : :
152 : : void
153 [ + + ]: 64500 : farbc( tk::Fields& U,
154 : : const std::vector< std::size_t >& farbcnodes,
155 : : const std::vector< tk::real >& farbcnorms )
156 : : // *****************************************************************************
157 : : // Set farfield boundary conditions at nodes
158 : : //! \param[in] U Solution vector at recent time step
159 : : //! \param[in] farbcnodes Nodes ids at which to set farfield BCs
160 : : //! \param[in] farbcnorms Normals at nodes at which to set farfield BCs
161 : : // *****************************************************************************
162 : : {
163 : : const auto& tf = g_cfg.get< tag::bc_far >();
164 [ + + ]: 64500 : if (tf.get< tag::sidesets >().empty()) return;
165 : :
166 : : Assert( farbcnodes.size()*3 == farbcnorms.size(), "Size mismatch" );
167 : :
168 : : // cppcheck-suppress unreadVariable
169 [ - + ]: 2425 : tk::real fr = tf.get< tag::density >();
170 : :
171 : : const auto& fue = tf.get< tag::velocity >();
172 [ - + ][ - - ]: 2425 : ErrChk( !fue.empty(), "No farfield velocity specified" );
[ - - ][ - - ]
173 : : // cppcheck-suppress unreadVariable
174 : 2425 : tk::real fu = fue[0];
175 : : // cppcheck-suppress unreadVariable
176 : 2425 : tk::real fv = fue[1];
177 : : // cppcheck-suppress unreadVariable
178 : 2425 : tk::real fw = fue[2];
179 : :
180 : 2425 : tk::real fp = tf.get< tag::pressure >();
181 : :
182 [ + + ]: 55840 : for (std::size_t i=0; i<farbcnodes.size(); ++i) {
183 : 53415 : auto p = farbcnodes[i];
184 : 53415 : auto nx = farbcnorms[i*3+0];
185 : 53415 : auto ny = farbcnorms[i*3+1];
186 [ - + ]: 53415 : auto nz = farbcnorms[i*3+2];
187 : : auto& r = U(p,0);
188 : : auto& ru = U(p,1);
189 : : auto& rv = U(p,2);
190 : : auto& rw = U(p,3);
191 : : auto& re = U(p,4);
192 : : //auto vn = (ru*nx + rv*ny + rw*nz)/r;
193 [ - + ]: 53415 : auto vn = fu*nx + fv*ny + fw*nz;
194 : : //auto a = eos::soundspeed( r,
195 : : // eos::pressure( re - 0.5*(ru*ru + rv*rv + rw*rw)/r ) );
196 : : auto a = eos::soundspeed( fr, fp );
197 : 53415 : auto M = vn / a;
198 [ - + ]: 53415 : if (M <= -1.0) {
199 : : // supersonic inflow, all characteristics from outside
200 : 0 : r = fr;
201 : 0 : ru = fr * fu;
202 : 0 : rv = fr * fv;
203 : 0 : rw = fr * fw;
204 : 0 : re = eos::totalenergy( fr, fu, fv, fw, fp );
205 [ + - ][ + + ]: 53415 : } else if (M > -1.0 && M < 0.0) {
206 : : // subsonic inflow: 1 outgoing and 4 incoming characteristics,
207 : : // pressure from inside, rest from outside
208 : 26697 : auto pr = eos::pressure( re - 0.5*(ru*ru + rv*rv + rw*rw)/r );
209 : 26697 : r = fr;
210 : 26697 : ru = fr * fu;
211 : 26697 : rv = fr * fv;
212 : 26697 : rw = fr * fw;
213 : 26697 : re = eos::totalenergy( fr, fu, fv, fw, pr );
214 [ + - ]: 53415 : } else if (M >= 0.0 && M < 1.0) {
215 : : // subsonic outflow: 1 incoming and 4 outgoing characteristics,
216 : : // pressure from outside, rest from inside
217 : 26718 : re = eos::totalenergy( r, ru/r, rv/r, rw/r, fp );
218 : : }
219 : : }
220 : : }
221 : :
222 : : void
223 : 64500 : prebc( tk::Fields& U,
224 : : const std::vector< std::size_t >& prebcnodes,
225 : : const std::vector< tk::real >& prebcvals )
226 : : // *****************************************************************************
227 : : // Set pressure boundary conditions at nodes
228 : : //! \param[in] U Solution vector at recent time step
229 : : //! \param[in] prebcnodes Node ids at which to set pressure BCs
230 : : //! \param[in] prebcvals Density and pressure values at pressure BC nodes
231 : : // *****************************************************************************
232 : : {
233 : : Assert( prebcnodes.size()*2 == prebcvals.size(), "Size mismatch" );
234 : :
235 [ + + ]: 74834 : for (std::size_t i=0; i<prebcnodes.size(); ++i) {
236 : 10334 : auto p = prebcnodes[i];
237 : 10334 : U(p,0) = prebcvals[i*2+0];
238 : 10334 : U(p,4) = eos::totalenergy( U(p,0), U(p,1)/U(p,0), U(p,2)/U(p,0),
239 : 10334 : U(p,3)/U(p,0), prebcvals[i*2+1] );
240 : : }
241 : 64500 : }
242 : :
243 : : } // physics::
|
