Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/Ghosts.cpp
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 Definitions file for generating ghost data structures
9 : : \details Definitions file for asynchronous distributed
10 : : ghost data structures using Charm++.
11 : : */
12 : : // *****************************************************************************
13 : :
14 : : #include "Ghosts.hpp"
15 : : #include "DerivedData.hpp"
16 : : #include "Reorder.hpp"
17 : : #include "Around.hpp"
18 : : #include "ChareStateCollector.hpp"
19 : :
20 : : extern tk::CProxy_ChareStateCollector stateProxy;
21 : :
22 : : using inciter::Ghosts;
23 : :
24 : 1099 : Ghosts::Ghosts( const CProxy_Discretization& disc,
25 : : const std::map< int, std::vector< std::size_t > >& bface,
26 : : const std::vector< std::size_t >& triinpoel,
27 : : std::size_t nunk,
28 : 1099 : CkCallback cbDone ) :
29 : : m_disc( disc ),
30 : : m_nunk( nunk ),
31 : : m_inpoel( Disc()->Inpoel() ),
32 : 1099 : m_coord( Disc()->Coord() ),
33 [ + - ]: 2198 : m_fd( m_inpoel, bface, tk::remap(triinpoel,Disc()->Lid()) ),
34 : 1099 : m_geoFace( tk::genGeoFaceTri( m_fd.Nipfac(), m_fd.Inpofa(), m_coord) ),
35 : 1099 : m_geoElem( tk::genGeoElemTet( m_inpoel, m_coord ) ),
36 : 1099 : m_nfac( m_fd.Inpofa().size()/3 ),
37 : : m_bndFace(),
38 : : m_sendGhost(),
39 : : m_ghost(),
40 : : m_exptGhost(),
41 : : m_bid(),
42 : : m_esup(),
43 : : m_initial( 1 ),
44 : : m_ncomfac( 0 ),
45 : : m_nadj( 0 ),
46 : : m_ncomEsup( 0 ),
47 : : m_ipface(),
48 : : m_ghostData(),
49 : : m_ghostReq( 0 ),
50 : : m_expChBndFace(),
51 : : m_infaces(),
52 : : m_esupc(),
53 [ + - ][ + - ]: 5495 : m_cbAfterDone( cbDone )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
54 : : // *****************************************************************************
55 : : // Constructor
56 : : //! \param[in] disc Discretization proxy
57 : : //! \param[in] bface Boundary-faces mapped to side set ids
58 : : //! \param[in] triinpoel Boundary-face connectivity
59 : : //! \param[in] nunk Number of unknowns
60 : : //! \param[in] cbDone Function to continue with when Ghosts have been computed
61 : : // *****************************************************************************
62 : : {
63 [ + + ][ + + ]: 2160 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
64 [ + + ]: 1061 : g_inputdeck.get< tag::cmd, tag::quiescence >())
65 [ + - ][ + - ]: 608 : stateProxy.ckLocalBranch()->insert( "Ghosts", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
[ + - ]
66 : : "Ghosts" );
67 : 1099 : }
68 : :
69 : : void
70 : 1099 : Ghosts::startCommSetup()
71 : : // *****************************************************************************
72 : : // Start setup of communication maps for cell-centered schemes
73 : : // *****************************************************************************
74 : : {
75 : : // Ensure that mesh partition is not leaky
76 [ - + ][ - - ]: 1099 : Assert( !tk::leakyPartition(m_fd.Esuel(), m_inpoel, m_coord),
[ - - ][ - - ]
77 : : "Input mesh to Ghosts leaky" );
78 : :
79 : : // Ensure mesh physical boundary for the entire problem not leaky,
80 : : // effectively checking if the user has specified boundary conditions on all
81 : : // physical boundary faces
82 : 1099 : bndIntegral();
83 : 1099 : }
84 : :
85 : : void
86 : 1099 : Ghosts::bndIntegral()
87 : : // *****************************************************************************
88 : : // Compute partial boundary surface integral and sum across all chares
89 : : //! \details This function computes a partial surface integral over the boundary
90 : : //! of the faces of this mesh partition then sends its contribution to perform
91 : : //! the integral acorss the total problem boundary. After the global sum a
92 : : //! non-zero vector result indicates a leak, e.g., a hole in the boundary
93 : : //! which indicates an error in the boundary face data structures used to
94 : : //! compute the partial surface integrals.
95 : : // *****************************************************************************
96 : : {
97 : : // Storage for surface integral over our mesh chunk physical boundary
98 [ + - ]: 1099 : std::vector< tk::real > s{{ 0.0, 0.0, 0.0 }};
99 : :
100 : : // Integrate over all physical boundary faces
101 [ + - ][ + + ]: 123789 : for (std::size_t f=0; f<m_fd.Nbfac(); ++f) {
102 [ + - ][ + - ]: 122690 : s[0] += m_geoFace(f,0) * m_geoFace(f,1);
103 [ + - ][ + - ]: 122690 : s[1] += m_geoFace(f,0) * m_geoFace(f,2);
104 [ + - ][ + - ]: 122690 : s[2] += m_geoFace(f,0) * m_geoFace(f,3);
105 : : }
106 : :
107 [ + - ]: 1099 : s.push_back( 1.0 ); // positive: call-back to resizeComm() after reduction
108 [ + - ][ + - ]: 1099 : s.push_back( static_cast< tk::real >( Disc()->MeshId() ) );
109 : :
110 : : // Send contribution to host summing partial surface integrals
111 [ + - ]: 1099 : contribute( s, CkReduction::sum_double,
112 [ + - ][ + - ]: 2198 : CkCallback(CkReductionTarget(Transporter,bndint), Disc()->Tr()) );
[ + - ]
113 : 1099 : }
114 : :
115 : : void
116 : 1099 : Ghosts::resizeComm()
117 : : // *****************************************************************************
118 : : // Start sizing communication buffers and setting up ghost data
119 : : // *****************************************************************************
120 : : {
121 : : // Enable SDAG wait for setting up chare boundary faces
122 [ + - ][ + - ]: 1099 : thisProxy[ thisIndex ].wait4fac();
123 : :
124 [ + - ]: 1099 : auto d = Disc();
125 : :
126 : 1099 : const auto& gid = d->Gid();
127 : 1099 : const auto& inpofa = m_fd.Inpofa();
128 : 1099 : const auto& esuel = m_fd.Esuel();
129 : :
130 : : // Perform leak test on mesh partition
131 [ + - ][ - + ]: 1099 : Assert( !tk::leakyPartition( esuel, m_inpoel, m_coord ),
[ - - ][ - - ]
[ - - ]
132 : : "Mesh partition leaky" );
133 : :
134 : : // Activate SDAG waits for face adjacency map (ghost data) calculation
135 [ + - ][ + - ]: 1099 : thisProxy[ thisIndex ].wait4ghost();
136 [ + - ][ + - ]: 1099 : thisProxy[ thisIndex ].wait4esup();
137 : :
138 : : // Invert inpofa to enable searching for faces based on (global) node triplets
139 [ - + ][ - - ]: 1099 : Assert( inpofa.size() % 3 == 0, "Inpofa must contain triplets" );
[ - - ][ - - ]
140 [ + + ]: 590328 : for (std::size_t f=0; f<inpofa.size()/3; ++f)
141 [ + - ]: 1178458 : m_ipface.insert( {{{ gid[ inpofa[f*3+0] ],
142 : 589229 : gid[ inpofa[f*3+1] ],
143 : 589229 : gid[ inpofa[f*3+2] ] }}} );
144 : :
145 : : // At this point ipface has node-id-triplets (faces) on the internal
146 : : // chare-domain and on the physical boundary but not on chare boundaries,
147 : : // hence the name internal + physical boundary faces.
148 : :
149 : : // Build a set of faces (each face given by 3 global node IDs) associated to
150 : : // chares we potentially share boundary faces with.
151 : 2198 : tk::UnsMesh::FaceSet potbndface;
152 [ + + ]: 274455 : for (std::size_t e=0; e<esuel.size()/4; ++e) { // for all our tets
153 : 273356 : auto mark = e*4;
154 [ + + ]: 1366780 : for (std::size_t f=0; f<4; ++f) // for all tet faces
155 [ + + ]: 1093424 : if (esuel[mark+f] == -1) { // if face has no outside-neighbor tet
156 : : // if does not exist among the internal and physical boundary faces,
157 : : // store as a potential chare-boundary face
158 : 160346 : tk::UnsMesh::Face t{{ gid[ m_inpoel[ mark + tk::lpofa[f][0] ] ],
159 : 160346 : gid[ m_inpoel[ mark + tk::lpofa[f][1] ] ],
160 : 320692 : gid[ m_inpoel[ mark + tk::lpofa[f][2] ] ] }};
161 [ + - ][ + + ]: 160346 : if (m_ipface.find(t) == end(m_ipface)) {
162 [ + - ][ - + ]: 37656 : Assert( m_expChBndFace.insert(t).second,
[ - - ][ - - ]
[ - - ]
163 : : "Store expected chare-boundary face" );
164 [ + - ]: 37656 : potbndface.insert( t );
165 : : }
166 : : }
167 : : }
168 : :
169 [ - + ][ - - ]: 1099 : if ( g_inputdeck.get< tag::cmd, tag::feedback >() ) d->Tr().chbndface();
170 : :
171 : : // In the following we assume that the size of the (potential) boundary-face
172 : : // adjacency map above does not necessarily equal to that of the node
173 : : // adjacency map. This is because while a node can be shared at a single
174 : : // corner or along an edge, that does not necessarily share a face as well
175 : : // (in other words, shared nodes or edges can exist that are not part of a
176 : : // shared face). So the chares we communicate with across faces are not
177 : : // necessarily the same as the chares we would communicate nodes with.
178 : : //
179 : : // Since the sizes of the node and face adjacency maps are not the same, while
180 : : // sending the faces on chare boundaries would be okay, however, the receiver
181 : : // would not necessarily know how many chares it must receive from. To solve
182 : : // this problem we send to chares which we share at least a single node with,
183 : : // i.e., rely on the node-adjacency map. Note that to all chares we share at
184 : : // least a single node with we send all our potential chare-boundary faces.
185 : : // This is the same list of faces to all chares we send.
186 : : //
187 : : // Another underlying assumption here is, of course, that the size of the face
188 : : // adjacency map is always smaller than or equal to that of the node adjacency
189 : : // map, which is always true. Since the receive side already knows how many
190 : : // fellow chares it must receive shared node ids from, we use that to detect
191 : : // completion of the number of receives in comfac(). This simplifies the
192 : : // communication pattern and code.
193 : :
194 : : // Send sets of faces adjacent to chare boundaries to fellow workers (if any)
195 [ + + ]: 1099 : if (d->NodeCommMap().empty()) // in serial, skip setting up ghosts altogether
196 [ + - ]: 38 : faceAdj();
197 : : else
198 : : // for all chares we share nodes with
199 [ + + ]: 11573 : for (const auto& c : d->NodeCommMap()) {
200 [ + - ][ + - ]: 10512 : thisProxy[ c.first ].comfac( thisIndex, potbndface );
201 : : }
202 : :
203 [ + - ]: 1099 : ownfac_complete();
204 : 1099 : }
205 : :
206 : : void
207 : 10512 : Ghosts::comfac( int fromch, const tk::UnsMesh::FaceSet& infaces )
208 : : // *****************************************************************************
209 : : // Receive unique set of faces we potentially share with/from another chare
210 : : //! \param[in] fromch Sender chare id
211 : : //! \param[in] infaces Unique set of faces we potentially share with fromch
212 : : // *****************************************************************************
213 : : {
214 [ + + ][ + + ]: 20558 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
215 [ + + ]: 10046 : g_inputdeck.get< tag::cmd, tag::quiescence >())
216 [ + - ][ + - ]: 5868 : stateProxy.ckLocalBranch()->insert( "Ghosts", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
217 : : "comfac" );
218 : :
219 : : // Buffer up incoming data
220 : 10512 : m_infaces[ fromch ] = infaces;
221 : :
222 : : // if we have heard from all fellow chares that we share at least a single
223 : : // node, edge, or face with
224 [ + + ]: 10512 : if (++m_ncomfac == Disc()->NodeCommMap().size()) {
225 : 1061 : m_ncomfac = 0;
226 : 1061 : comfac_complete();
227 : : }
228 : 10512 : }
229 : :
230 : : void
231 : 1061 : Ghosts::bndFaces()
232 : : // *****************************************************************************
233 : : // Compute chare-boundary faces
234 : : //! \details This is called when both send and receives are completed on a
235 : : //! chare and thus we are ready to compute chare-boundary faces and ghost data.
236 : : // *****************************************************************************
237 : : {
238 : 1061 : auto d = Disc();
239 [ - + ]: 1061 : if ( g_inputdeck.get< tag::cmd, tag::feedback >() ) d->Tr().chcomfac();
240 : 1061 : const auto& esuel = m_fd.Esuel();
241 : 1061 : const auto& gid = d->Gid();
242 : :
243 [ + + ]: 11573 : for (const auto& in : m_infaces) {
244 : : // Find sender chare among chares we potentially share faces with. Note that
245 : : // it is feasible that a sender chare called us but we do not have a set of
246 : : // faces associated to that chare. This can happen if we only share a single
247 : : // node or an edge but not a face with that chare.
248 [ + - ]: 10512 : auto& bndface = m_bndFace[ in.first ]; // will associate to sender chare
249 : : // Try to find incoming faces on our chare boundary with other chares. If
250 : : // found, generate and assign new local face ID, associated to sender chare.
251 [ + + ]: 656879 : for (std::size_t e=0; e<esuel.size()/4; ++e) { // for all our tets
252 : 646367 : auto mark = e*4;
253 [ + + ]: 3231835 : for (std::size_t f=0; f<4; ++f) { // for all cell faces
254 [ + + ]: 2585468 : if (esuel[mark+f] == -1) { // if face has no outside-neighbor tet
255 : 635138 : tk::UnsMesh::Face t{{ gid[ m_inpoel[ mark + tk::lpofa[f][0] ] ],
256 : 635138 : gid[ m_inpoel[ mark + tk::lpofa[f][1] ] ],
257 : 1270276 : gid[ m_inpoel[ mark + tk::lpofa[f][2] ] ] }};
258 : : // if found among the incoming faces and if not one of our internal
259 : : // nor physical boundary faces
260 [ + - ][ + + ]: 672794 : if ( in.second.find(t) != end(in.second) &&
[ + - ]
261 [ + - ][ + + ]: 672794 : m_ipface.find(t) == end(m_ipface) ) {
262 [ + - ]: 37656 : bndface[t][0] = m_nfac++; // assign new local face ID
263 : : }
264 : : }
265 : : }
266 : : }
267 : : // If at this point if we have not found any face among our faces we
268 : : // potentially share with fromch, there is no need to keep an empty set of
269 : : // faces associated to fromch as we only share nodes or edges with it, but
270 : : // not faces.
271 [ + + ][ + - ]: 10512 : if (bndface.empty()) m_bndFace.erase( in.first );
272 : : }
273 : :
274 : 1061 : tk::destroy(m_ipface);
275 : 1061 : tk::destroy(m_infaces);
276 : :
277 : : // Ensure all expected faces have been received
278 [ - + ][ - - ]: 1061 : Assert( receivedChBndFaces(),
[ - - ][ - - ]
279 : : "Expected and received chare boundary faces mismatch" );
280 : :
281 : : // Basic error checking on chare-boundary-face map
282 [ + - ][ - + ]: 1061 : Assert( m_bndFace.find( thisIndex ) == m_bndFace.cend(),
[ - - ][ - - ]
[ - - ]
283 : : "Face-communication map should not contain data for own chare ID" );
284 : :
285 : : // Store (local) tet ID adjacent to our chare boundary from the inside
286 [ + + ]: 147079 : for (std::size_t e=0; e<esuel.size()/4; ++e) { // for all our tets
287 : 146018 : auto mark = e*4;
288 [ + + ]: 730090 : for (std::size_t f=0; f<4; ++f) { // for all cell faces
289 [ + + ]: 584072 : if (esuel[mark+f] == -1) { // if face has no outside-neighbor tet
290 : 107840 : tk::UnsMesh::Face t{{ gid[ m_inpoel[ mark + tk::lpofa[f][0] ] ],
291 : 107840 : gid[ m_inpoel[ mark + tk::lpofa[f][1] ] ],
292 : 215680 : gid[ m_inpoel[ mark + tk::lpofa[f][2] ] ] }};
293 [ + - ]: 107840 : auto c = findchare(t);
294 [ + + ]: 107840 : if (c > -1) {
295 [ + - ]: 37656 : auto& lbndface = tk::ref_find( m_bndFace, c );
296 [ + - ]: 37656 : auto& face = tk::ref_find( lbndface, t );
297 : 37656 : face[1] = e; // store (local) inner tet ID adjacent to face
298 : : }
299 : : }
300 : : }
301 : : }
302 : :
303 : : // At this point m_bndFace is complete on this PE. This means that starting
304 : : // from the sets of faces we potentially share with fellow chares we now
305 : : // only have those faces we actually share faces with (through which we need
306 : : // to communicate later). Also, m_bndFace not only has the unique faces
307 : : // associated to fellow chares, but also a newly assigned local face ID as
308 : : // well as the local id of the inner tet adjacent to the face. Continue by
309 : : // starting setting up ghost data
310 : 1061 : setupGhost();
311 : : // Besides setting up our own ghost data, we also issue requests (for ghost
312 : : // data) to those chares which we share faces with. Note that similar to
313 : : // comfac() we are calling reqGhost() by going through the node communication
314 : : // map instead, which may send requests to those chare we do not share faces
315 : : // with. This is so that we can test for completing by querying the size of
316 : : // the already complete node commincation map in reqGhost. Requests in
317 : : // sendGhost will only be fullfilled based on m_ghostData.
318 [ + + ]: 11573 : for (const auto& c : d->NodeCommMap()) // for all chares we share nodes with
319 [ + - ][ + - ]: 10512 : thisProxy[ c.first ].reqGhost();
320 : 1061 : }
321 : :
322 : : void
323 : 1061 : Ghosts::setupGhost()
324 : : // *****************************************************************************
325 : : // Setup own ghost data on this chare
326 : : // *****************************************************************************
327 : : {
328 : 1061 : auto d = Disc();
329 : 1061 : const auto& gid = d->Gid();
330 : :
331 : : // Enlarge elements surrounding faces data structure for ghosts
332 [ + - ]: 1061 : m_fd.Esuf().resize( 2*m_nfac, -2 );
333 [ + - ]: 1061 : m_fd.Inpofa().resize( 3*m_nfac, 0 );
334 : : // Enlarge face geometry data structure for ghosts
335 : 1061 : m_geoFace.resize( m_nfac, 0.0 );
336 : :
337 : 1061 : const auto& esuel = m_fd.Esuel();
338 : :
339 : : // Collect tet ids, their face connectivity (given by 3 global node IDs, each
340 : : // triplet for potentially multiple faces on the chare boundary), and their
341 : : // elem geometry data (see GhostData) associated to fellow chares adjacent to
342 : : // chare boundaries. Once received by fellow chares, these tets will become
343 : : // known as ghost elements and their data as ghost data.
344 [ + + ]: 147079 : for (std::size_t e=0; e<esuel.size()/4; ++e) { // for all our tets
345 : 146018 : auto mark = e*4;
346 [ + + ]: 730090 : for (std::size_t f=0; f<4; ++f) { // for all cell faces
347 [ + + ]: 584072 : if (esuel[mark+f] == -1) { // if face has no outside-neighbor tet
348 : 107840 : tk::UnsMesh::Face t{{ gid[ m_inpoel[ mark + tk::lpofa[f][0] ] ],
349 : 107840 : gid[ m_inpoel[ mark + tk::lpofa[f][1] ] ],
350 : 215680 : gid[ m_inpoel[ mark + tk::lpofa[f][2] ] ] }};
351 [ + - ]: 107840 : auto c = findchare(t);
352 : : // It is possible that we do not find the chare for this face. We are
353 : : // looping through all of our tets and interrogating all faces that do
354 : : // not have neighboring tets but we only care about chare-boundary faces
355 : : // here as only those need ghost data. (esuel may also contain
356 : : // physical boundary faces)
357 [ + + ]: 107840 : if (c > -1) {
358 : : // Will store ghost data associated to neighbor chare
359 [ + - ]: 37656 : auto& ghost = m_ghostData[ c ];
360 : : // Store tet id adjacent to chare boundary as key for ghost data
361 [ + - ]: 37656 : auto& tuple = ghost[ e ];
362 : : // If tetid e has not yet been encountered, store geometry (only once)
363 : 37656 : auto& nodes = std::get< 0 >( tuple );
364 [ + + ]: 37656 : if (nodes.empty()) {
365 [ + - ]: 29938 : std::get< 1 >( tuple ) = m_geoElem[ e ];
366 : :
367 : 29938 : auto& ncoord = std::get< 2 >( tuple );
368 : 29938 : ncoord[0] = m_coord[0][ m_inpoel[ mark+f ] ];
369 : 29938 : ncoord[1] = m_coord[1][ m_inpoel[ mark+f ] ];
370 : 29938 : ncoord[2] = m_coord[2][ m_inpoel[ mark+f ] ];
371 : :
372 : 29938 : std::get< 3 >( tuple ) = f;
373 : :
374 : 59876 : std::get< 4 >( tuple ) = {{ gid[ m_inpoel[ mark ] ],
375 : 29938 : gid[ m_inpoel[ mark+1 ] ],
376 : 29938 : gid[ m_inpoel[ mark+2 ] ],
377 : 29938 : gid[ m_inpoel[ mark+3 ] ] }};
378 : : }
379 : : // (Always) store face node IDs on chare boundary, even if tetid e has
380 : : // already been stored. Thus we store potentially multiple faces along
381 : : // the same chare-boundary. This happens, e.g., when the boundary
382 : : // between chares is zig-zaggy enough to have 2 or even 3 faces of the
383 : : // same tet.
384 [ + - ]: 37656 : nodes.push_back( t[0] );
385 [ + - ]: 37656 : nodes.push_back( t[1] );
386 [ + - ]: 37656 : nodes.push_back( t[2] );
387 [ - + ][ - - ]: 37656 : Assert( nodes.size() <= 4*3, "Overflow of faces/tet to send" );
[ - - ][ - - ]
388 : : }
389 : : }
390 : : }
391 : : }
392 : :
393 : : // Basic error checking on local ghost data
394 [ + - ][ - + ]: 1061 : Assert( m_ghostData.find( thisIndex ) == m_ghostData.cend(),
[ - - ][ - - ]
[ - - ]
395 : : "Chare-node adjacency map should not contain data for own chare ID" );
396 : :
397 : : // More in-depth error checking on local ghost data
398 [ + + ]: 6387 : for (const auto& c : m_ghostData)
399 [ + + ]: 35264 : for ([[maybe_unused]] const auto& t : c.second) {
400 [ - + ][ - - ]: 29938 : Assert( !std::get< 0 >( t.second ).empty(),
[ - - ][ - - ]
401 : : "Emtpy face vector in ghost data" );
402 [ - + ][ - - ]: 29938 : Assert( std::get< 0 >( t.second ).size() % 3 == 0,
[ - - ][ - - ]
403 : : "Face node IDs must be triplets" );
404 [ - + ][ - - ]: 29938 : Assert( std::get< 0 >( t.second ).size() <= 4*3, // <= 4*3 (4*numfaces)
[ - - ][ - - ]
405 : : "Max number of faces for a single ghost tet is 4" );
406 [ - + ][ - - ]: 29938 : Assert( !std::get< 1 >( t.second ).empty(),
[ - - ][ - - ]
407 : : "No elem geometry data for ghost" );
408 [ - + ][ - - ]: 29938 : Assert( std::get< 1 >( t.second ).size() == m_geoElem.nprop(),
[ - - ][ - - ]
409 : : "Elem geometry data for ghost must be for single tet" );
410 [ - + ][ - - ]: 29938 : Assert( !std::get< 2 >( t.second ).empty(),
[ - - ][ - - ]
411 : : "No nodal coordinate data for ghost" );
412 : : }
413 : :
414 : 1061 : ownghost_complete();
415 : 1061 : }
416 : :
417 : : void
418 : 10512 : Ghosts::reqGhost()
419 : : // *****************************************************************************
420 : : // Receive requests for ghost data
421 : : // *****************************************************************************
422 : : {
423 [ + + ][ + + ]: 20558 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
424 [ + + ]: 10046 : g_inputdeck.get< tag::cmd, tag::quiescence >())
425 [ + - ][ + - ]: 5868 : stateProxy.ckLocalBranch()->insert( "Ghosts", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
426 : : "reqGhost" );
427 : :
428 : : // If every chare we communicate with has requested ghost data from us, we may
429 : : // fulfill the requests, but only if we have already setup our ghost data.
430 [ + + ]: 10512 : if (++m_ghostReq == Disc()->NodeCommMap().size()) {
431 : 1061 : m_ghostReq = 0;
432 : 1061 : reqghost_complete();
433 : : }
434 : 10512 : }
435 : :
436 : : void
437 : 1061 : Ghosts::sendGhost()
438 : : // *****************************************************************************
439 : : // Send all of our ghost data to fellow chares
440 : : // *****************************************************************************
441 : : {
442 [ + + ][ + + ]: 2084 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
443 [ + + ]: 1023 : g_inputdeck.get< tag::cmd, tag::quiescence >())
444 [ + - ][ + - ]: 572 : stateProxy.ckLocalBranch()->insert( "Ghosts", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
445 : : "sendGhost" );
446 : :
447 [ + + ]: 6387 : for (const auto& c : m_ghostData)
448 [ + - ][ + - ]: 5326 : thisProxy[ c.first ].comGhost( thisIndex, c.second );
449 : :
450 [ - + ]: 1061 : if ( g_inputdeck.get< tag::cmd, tag::feedback >() ) Disc()->Tr().chghost();
451 : 1061 : }
452 : :
453 : : void
454 : 5326 : Ghosts::comGhost( int fromch, const GhostData& ghost )
455 : : // *****************************************************************************
456 : : // Receive ghost data on chare boundaries from fellow chare
457 : : //! \param[in] fromch Caller chare ID
458 : : //! \param[in] ghost Ghost data, see Inciter/FaceData.h for the type
459 : : // *****************************************************************************
460 : : {
461 [ + + ][ + + ]: 10412 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
462 [ + + ]: 5086 : g_inputdeck.get< tag::cmd, tag::quiescence >())
463 [ + - ][ + - ]: 2862 : stateProxy.ckLocalBranch()->insert( "Ghosts", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
464 : : "comGhost" );
465 : :
466 : 5326 : auto d = Disc();
467 : 5326 : const auto& lid = d->Lid();
468 : 5326 : auto& inpofa = m_fd.Inpofa();
469 : 5326 : auto ncoord = m_coord[0].size();
470 : :
471 : : // nodelist with fromch, currently only used for an assert
472 : 5326 : [[maybe_unused]] const auto& nl = tk::cref_find( d->NodeCommMap(), fromch );
473 : :
474 : 5326 : auto& ghostelem = m_ghost[ fromch ]; // will associate to sender chare
475 : :
476 : : // Store ghost data coming from chare
477 [ + + ]: 35264 : for (const auto& g : ghost) { // loop over incoming ghost data
478 : 29938 : auto e = g.first; // remote/ghost tet id outside of chare boundary
479 : 29938 : const auto& nodes = std::get< 0 >( g.second ); // node IDs of face(s)
480 : 29938 : const auto& geo = std::get< 1 >( g.second ); // ghost elem geometry data
481 : 29938 : const auto& coordg = std::get< 2 >( g.second ); // coordinate of ghost node
482 : 29938 : const auto& inpoelg = std::get< 4 >( g.second ); // inpoel of ghost tet
483 : :
484 [ - + ][ - - ]: 29938 : Assert( nodes.size() % 3 == 0, "Face node IDs must be triplets" );
[ - - ][ - - ]
485 [ - + ][ - - ]: 29938 : Assert( nodes.size() <= 4*3, "Overflow of faces/tet received" );
[ - - ][ - - ]
486 [ - + ][ - - ]: 29938 : Assert( geo.size() % 5 == 0, "Ghost geometry size mismatch" );
[ - - ][ - - ]
487 [ - + ][ - - ]: 29938 : Assert( geo.size() == m_geoElem.nprop(), "Ghost geometry number mismatch" );
[ - - ][ - - ]
488 [ - + ][ - - ]: 29938 : Assert( coordg.size() == 3, "Incorrect ghost node coordinate size" );
[ - - ][ - - ]
489 [ - + ][ - - ]: 29938 : Assert( inpoelg.size() == 4, "Incorrect ghost inpoel size" );
[ - - ][ - - ]
490 : :
491 [ + + ]: 67594 : for (std::size_t n=0; n<nodes.size()/3; ++n) { // face(s) of ghost e
492 : : // node IDs of face on chare boundary
493 : 37656 : tk::UnsMesh::Face t{{ nodes[n*3+0], nodes[n*3+1], nodes[n*3+2] }};
494 : : // must find t in nodelist of chare-boundary adjacent to fromch
495 [ + - ][ + - ]: 37656 : Assert( nl.find(t[0]) != end(nl) &&
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
[ - - ]
496 : : nl.find(t[1]) != end(nl) &&
497 : : nl.find(t[2]) != end(nl),
498 : : "Ghost face not found in chare-node adjacency map on receiving end" );
499 : : // must find face in boundary-face adjacency map for fromch
500 [ + - ][ + - ]: 37656 : Assert( tk::cref_find(m_bndFace,fromch).find( t ) !=
[ + - ][ - + ]
[ - - ][ - - ]
[ - - ]
501 : : tk::cref_find(m_bndFace,fromch).cend(), "Ghost face not "
502 : : "found in boundary-face adjacency map on receiving end" );
503 : : // find local face & tet ids for t
504 [ + - ][ + - ]: 37656 : auto id = tk::cref_find( tk::cref_find(m_bndFace,fromch), t );
505 : : // compute face geometry for chare-boundary face
506 [ + - ]: 37656 : addGeoFace(t, id);
507 : : // add node-triplet to node-face connectivity
508 [ + - ]: 37656 : inpofa[3*id[0]+0] = tk::cref_find( lid, t[2] );
509 [ + - ]: 37656 : inpofa[3*id[0]+1] = tk::cref_find( lid, t[1] );
510 [ + - ]: 37656 : inpofa[3*id[0]+2] = tk::cref_find( lid, t[0] );
511 : :
512 : : // if ghost tet id not yet encountered on boundary with fromch
513 [ + - ]: 37656 : auto i = ghostelem.find( e );
514 [ + + ]: 37656 : if (i != end(ghostelem)) {
515 : : // fill in elements surrounding face
516 [ + - ]: 7718 : addEsuf(id, i->second);
517 : : // fill in elements surrounding element
518 [ + - ]: 7718 : addEsuel(id, i->second, t);
519 : : } else {
520 : : // fill in elements surrounding face
521 [ + - ]: 29938 : addEsuf(id, m_nunk);
522 : : // fill in elements surrounding element
523 [ + - ]: 29938 : addEsuel(id, m_nunk, t);
524 [ + - ]: 29938 : ghostelem[e] = m_nunk; // assign new local tet id to remote ghost id
525 [ + - ]: 29938 : m_geoElem.push_back( geo );// store ghost elem geometry
526 : 29938 : ++m_nunk; // increase number of unknowns on this chare
527 : 29938 : std::size_t counter = 0;
528 [ + + ]: 149690 : for (std::size_t gp=0; gp<4; ++gp) {
529 [ + - ]: 119752 : auto it = lid.find( inpoelg[gp] );
530 : : std::size_t lp;
531 [ + + ]: 119752 : if (it != end(lid))
532 : 99189 : lp = it->second;
533 : : else {
534 [ - + ][ - - ]: 20563 : Assert( nodes.size() == 3, "Expected node not found in lid" );
[ - - ][ - - ]
535 [ - + ][ - - ]: 20563 : Assert( gp == std::get< 3 >( g.second ),
[ - - ][ - - ]
536 : : "Ghost node not matching correct entry in ghost inpoel" );
537 : 20563 : lp = ncoord;
538 : 20563 : ++counter;
539 : : }
540 [ + - ]: 119752 : m_inpoel.push_back( lp ); // store ghost element connectivity
541 : : }
542 : : // only a single or no ghost node should be found
543 [ - + ][ - - ]: 29938 : Assert( counter <= 1, "Incorrect number of ghost nodes detected. "
[ - - ][ - - ]
[ - - ][ - - ]
544 : : "Detected "+ std::to_string(counter) +" ghost nodes" );
545 [ + + ]: 29938 : if (counter == 1) {
546 [ + - ]: 20563 : m_coord[0].push_back( coordg[0] ); // store ghost node coordinate
547 [ + - ]: 20563 : m_coord[1].push_back( coordg[1] );
548 [ + - ]: 20563 : m_coord[2].push_back( coordg[2] );
549 [ - + ][ - - ]: 20563 : Assert( m_inpoel[ 4*(m_nunk-1)+std::get< 3 >( g.second ) ] == ncoord,
[ - - ][ - - ]
550 : : "Mismatch in extended inpoel for ghost element" );
551 : 20563 : ++ncoord; // increase number of nodes on this chare
552 : : }
553 : : }
554 : :
555 : : // additional tests to ensure that entries in inpoel and t/inpofa match
556 [ + - ][ - + ]: 37656 : Assert( nodetripletMatch(id, t) == 3,
[ - - ][ - - ]
[ - - ]
557 : : "Mismatch/Overmatch in inpoel and inpofa at chare-boundary face" );
558 : : }
559 : : }
560 : :
561 : : // Signal the runtime system that all workers have received their
562 : : // face-adjacency
563 [ + + ]: 5326 : if (++m_nadj == m_ghostData.size()) faceAdj();
564 : 5326 : }
565 : :
566 : : void
567 : 1099 : Ghosts::faceAdj()
568 : : // *****************************************************************************
569 : : // Continue after face adjacency communication map completed on this chare
570 : : //! \details At this point the face communication map has been established
571 : : //! on this chare. Proceed to set up the nodal-comm map.
572 : : // *****************************************************************************
573 : : {
574 : 1099 : m_nadj = 0;
575 : :
576 : 1099 : tk::destroy(m_bndFace);
577 : :
578 : : // Ensure that all elements surrounding faces (are correct) including those at
579 : : // chare boundaries
580 [ + + ]: 627984 : for (std::size_t f=0; f<m_nfac; ++f) {
581 [ - + ][ - - ]: 626885 : Assert( m_fd.Esuf()[2*f] > -1,
[ - - ][ - - ]
582 : : "Left element in esuf cannot be physical ghost" );
583 [ + - ][ + + ]: 626885 : if (f >= m_fd.Nbfac())
584 [ - + ][ - - ]: 504195 : Assert( m_fd.Esuf()[2*f+1] > -1,
[ - - ][ - - ]
585 : : "Right element in esuf for internal/chare faces cannot be a ghost" );
586 : : }
587 : :
588 : : // Ensure that all elements surrounding elements are correct including those
589 : : // at chare boundaries
590 : 1099 : const auto& esuel = m_fd.Esuel();
591 : 1099 : std::size_t nbound = 0;
592 [ + + ]: 274455 : for (std::size_t e=0; e<esuel.size()/4; ++e) {
593 [ + + ]: 1366780 : for (std::size_t f=0; f<4; ++f)
594 [ + + ]: 1093424 : if (esuel[4*e+f] == -1) ++nbound;
595 : : }
596 [ + - ][ - + ]: 1099 : Assert( nbound == m_fd.Nbfac(), "Incorrect number of ghost-element -1's in "
[ - - ][ - - ]
[ - - ]
597 : : "updated esuel" );
598 : :
599 : : // Error checking on ghost data
600 [ + + ]: 6425 : for(const auto& n : m_ghostData)
601 [ + + ]: 35264 : for([[maybe_unused]] const auto& i : n.second)
602 [ - + ][ - - ]: 29938 : Assert( i.first < m_fd.Esuel().size()/4, "Sender contains ghost tet id " );
[ - - ][ - - ]
603 : :
604 : : // Perform leak test on face geometry data structure enlarged by ghosts
605 [ + - ][ - + ]: 1099 : Assert( !leakyAdjacency(), "Face adjacency leaky" );
[ - - ][ - - ]
[ - - ]
606 [ + - ][ - + ]: 1099 : Assert( faceMatch(), "Chare-boundary element-face "
[ - - ][ - - ]
[ - - ]
607 : : "connectivity (esuf) does not match" );
608 : :
609 : : // Create new map of elements along chare boundary which are ghosts for
610 : : // neighboring chare, associated with that chare ID
611 [ + + ]: 6425 : for (const auto& [cid, cgd] : m_ghostData)
612 : : {
613 [ + - ]: 5326 : auto& sg = m_sendGhost[cid];
614 [ + + ]: 35264 : for (const auto& e : cgd)
615 : : {
616 [ + - ][ - + ]: 29938 : Assert(sg.find(e.first) == sg.end(), "Repeating element found in "
[ - - ][ - - ]
[ - - ]
617 : : "ghost data");
618 [ + - ]: 29938 : sg.insert(e.first);
619 : : }
620 [ - + ][ - - ]: 5326 : Assert(sg.size() == cgd.size(), "Incorrect size for sendGhost");
[ - - ][ - - ]
621 : : }
622 [ - + ][ - - ]: 1099 : Assert(m_sendGhost.size() == m_ghostData.size(), "Incorrect number of "
[ - - ][ - - ]
623 : : "chares in sendGhost");
624 : :
625 : : // Error checking on ghost data
626 [ + + ]: 6425 : for(const auto& n : m_sendGhost)
627 [ + + ]: 35264 : for([[maybe_unused]] const auto& i : n.second)
628 [ - + ][ - - ]: 29938 : Assert( i < m_fd.Esuel().size()/4, "Sender contains ghost tet id. " );
[ - - ][ - - ]
629 : :
630 : : // Generate and store Esup data-structure in a map
631 [ + - ]: 1099 : auto esup = tk::genEsup(m_inpoel, 4);
632 [ + - ][ + + ]: 96660 : for (std::size_t p=0; p<Disc()->Gid().size(); ++p)
633 : : {
634 [ + + ]: 1288174 : for (auto e : tk::Around(esup, p))
635 : : {
636 : : // since inpoel has been augmented with the face-ghost cell previously,
637 : : // esup also contains cells which are not on this mesh-chunk, hence the
638 : : // following test
639 [ + + ][ + - ]: 1192613 : if (e < m_fd.Esuel().size()/4) m_esup[p].push_back(e);
[ + - ]
640 : : }
641 : : }
642 : :
643 : : // Error checking on Esup map
644 [ + + ]: 96660 : for(const auto& p : m_esup)
645 [ + + ]: 1188985 : for([[maybe_unused]] const auto& e : p.second)
646 [ - + ][ - - ]: 1093424 : Assert( e < m_fd.Esuel().size()/4, "Esup contains tet id greater than "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
647 : : + std::to_string(m_fd.Esuel().size()/4-1) +" : "+ std::to_string(e) );
648 : :
649 [ + - ]: 1099 : auto meshid = Disc()->MeshId();
650 [ + - ]: 1099 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
651 [ + - ][ + - ]: 2198 : CkCallback(CkReductionTarget(Transporter,startEsup), Disc()->Tr()) );
[ + - ]
652 : 1099 : }
653 : :
654 : : void
655 : 1099 : Ghosts::nodeNeighSetup()
656 : : // *****************************************************************************
657 : : // Setup node-neighborhood (esup)
658 : : //! \details At this point the face-ghost communication map has been established
659 : : //! on this chare. This function begins generating the node-ghost comm map.
660 : : // *****************************************************************************
661 : : {
662 [ + + ]: 1099 : if (Disc()->NodeCommMap().empty())
663 : : // in serial, skip setting up node-neighborhood
664 : 38 : { comesup_complete(); }
665 : : else
666 : : {
667 : 1061 : const auto& nodeCommMap = Disc()->NodeCommMap();
668 : :
669 : : // send out node-neighborhood map
670 [ + + ]: 11573 : for (const auto& [cid, nlist] : nodeCommMap)
671 : : {
672 : 21024 : std::unordered_map< std::size_t, std::vector< std::size_t > > bndEsup;
673 : 10512 : std::unordered_map< std::size_t, std::vector< tk::real > > nodeBndCells;
674 [ + + ]: 69836 : for (const auto& p : nlist)
675 : : {
676 [ + - ][ + - ]: 59324 : auto pl = tk::cref_find(Disc()->Lid(), p);
677 : : // fill in the esup for the chare-boundary
678 [ + - ]: 59324 : const auto& pesup = tk::cref_find(m_esup, pl);
679 [ + - ][ + - ]: 59324 : bndEsup[p] = pesup;
680 : :
681 : : // fill a map with the element ids from esup as keys and geoElem as
682 : : // values, and another map containing these elements associated with
683 : : // the chare id with which they are node-neighbors.
684 [ + + ]: 360118 : for (const auto& e : pesup)
685 : : {
686 [ + - ][ + - ]: 300794 : nodeBndCells[e] = m_geoElem[e];
687 : :
688 : : // add these esup-elements into map of elements along chare boundary
689 [ - + ][ - - ]: 300794 : Assert( e < m_fd.Esuel().size()/4, "Sender contains ghost tet id." );
[ - - ][ - - ]
690 [ + - ][ + - ]: 300794 : m_sendGhost[cid].insert(e);
691 : : }
692 : : }
693 : :
694 [ + - ][ + - ]: 10512 : thisProxy[cid].comEsup(thisIndex, bndEsup, nodeBndCells);
695 : : }
696 : : }
697 : :
698 : 1099 : ownesup_complete();
699 : 1099 : }
700 : :
701 : : void
702 : 10512 : Ghosts::comEsup( int fromch,
703 : : const std::unordered_map< std::size_t, std::vector< std::size_t > >& bndEsup,
704 : : const std::unordered_map< std::size_t, std::vector< tk::real > >&
705 : : nodeBndCells )
706 : : // *****************************************************************************
707 : : //! \brief Receive elements-surrounding-points data-structure for points on
708 : : // common boundary between receiving and sending neighbor chare, and the
709 : : // element geometries for these new elements
710 : : //! \param[in] fromch Sender chare id
711 : : //! \param[in] bndEsup Elements-surrounding-points data-structure from fromch
712 : : //! \param[in] nodeBndCells Map containing element geometries associated with
713 : : //! remote element IDs in the esup
714 : : // *****************************************************************************
715 : : {
716 : 10512 : auto& chghost = m_ghost[fromch];
717 : :
718 : : // Extend remote-local element id map and element geometry array
719 [ + + ]: 169441 : for (const auto& e : nodeBndCells)
720 : : {
721 : : // need to check following, because 'e' could have been added previously in
722 : : // remote-local element id map as a part of face-communication, i.e. as a
723 : : // face-ghost element
724 [ + - ][ + + ]: 158929 : if (chghost.find(e.first) == chghost.end())
725 : : {
726 [ + - ]: 128991 : chghost[e.first] = m_nunk;
727 [ + - ]: 128991 : m_geoElem.push_back(e.second);
728 : 128991 : ++m_nunk;
729 : : }
730 : : }
731 : :
732 : : // Store incoming data in comm-map buffer for Esup
733 [ + + ]: 69836 : for (const auto& [node, elist] : bndEsup)
734 : : {
735 [ + - ][ + - ]: 59324 : auto pl = tk::cref_find(Disc()->Lid(), node);
736 [ + - ]: 59324 : auto& pesup = m_esupc[pl];
737 [ + + ]: 360118 : for (auto e : elist)
738 : : {
739 [ + - ]: 300794 : auto el = tk::cref_find(chghost, e);
740 [ + - ]: 300794 : pesup.push_back(el);
741 : : }
742 : : }
743 : :
744 : : // if we have heard from all fellow chares that we share at least a single
745 : : // node, edge, or face with
746 [ + + ]: 10512 : if (++m_ncomEsup == Disc()->NodeCommMap().size()) {
747 : 1061 : m_ncomEsup = 0;
748 : 1061 : comesup_complete();
749 : : }
750 : 10512 : }
751 : :
752 : : void
753 : 1099 : Ghosts::adj()
754 : : // *****************************************************************************
755 : : // Finish up with adjacency maps, and do a global-sync to begin problem setup
756 : : //! \details At this point, the nodal- and face-adjacency has been set up. This
757 : : // function does some error checking on the nodal-adjacency and prepares
758 : : // for problem setup.
759 : : // *****************************************************************************
760 : : {
761 : : // combine own and communicated contributions to elements surrounding points
762 [ + + ]: 31608 : for (auto& [p, elist] : m_esupc)
763 : : {
764 [ + - ]: 30509 : auto& pesup = tk::ref_find(m_esup, p);
765 [ + + ]: 331303 : for ([[maybe_unused]] auto e : elist)
766 : : {
767 [ - + ][ - - ]: 300794 : Assert( e >= m_fd.Esuel().size()/4, "Non-ghost element received from "
[ - - ][ - - ]
768 : : "esup buffer." );
769 : : }
770 [ + - ]: 30509 : tk::concat< std::size_t >(std::move(elist), pesup);
771 : : }
772 : :
773 : 1099 : tk::destroy(m_ghostData);
774 : 1099 : tk::destroy(m_esupc);
775 : :
776 [ - + ][ - - ]: 1099 : if ( g_inputdeck.get< tag::cmd, tag::feedback >() ) Disc()->Tr().chadj();
[ - - ]
777 : :
778 : : // Error checking on ghost data
779 [ + + ]: 11611 : for(const auto& n : m_sendGhost)
780 [ + + ]: 169441 : for([[maybe_unused]] const auto& i : n.second)
781 [ - + ][ - - ]: 158929 : Assert( i < m_fd.Esuel().size()/4, "Sender contains ghost tet id. ");
[ - - ][ - - ]
782 : :
783 : : // Create a mapping between local ghost tet ids and zero-based boundary ids
784 [ + - ]: 2198 : std::vector< std::size_t > c( tk::sumvalsize( m_ghost ) );
785 : 1099 : std::size_t j = 0;
786 [ + + ]: 11611 : for (const auto& n : m_ghost) {
787 [ + + ]: 169441 : for(const auto& i : n.second) {
788 : 158929 : c[j++] = i.second;
789 : : }
790 : : }
791 [ + - ]: 1099 : m_bid = tk::assignLid( c );
792 : :
793 : : // Basic error checking on ghost tet ID map
794 [ + - ][ - + ]: 1099 : Assert( m_ghost.find( thisIndex ) == m_ghost.cend(),
[ - - ][ - - ]
[ - - ]
795 : : "Ghost id map should not contain data for own chare ID" );
796 : :
797 : : // Store expected ghost tet IDs
798 [ + + ]: 11611 : for (const auto& n : m_ghost)
799 [ + + ]: 169441 : for ([[maybe_unused]] const auto& g : n.second)
800 [ + - ][ - + ]: 158929 : Assert( m_exptGhost.insert( g.second ).second,
[ - - ][ - - ]
[ - - ]
801 : : "Failed to store local tetid as exptected ghost id" );
802 : :
803 : : // Callback function from DG/FV after ghost-setup is done
804 [ + - ]: 1099 : m_cbAfterDone.send();
805 : 1099 : }
806 : :
807 : : bool
808 : 1099 : Ghosts::leakyAdjacency()
809 : : // *****************************************************************************
810 : : // Perform leak-test on chare boundary faces
811 : : //! \details This function computes a surface integral over the boundary of the
812 : : //! faces after the face adjacency communication map is completed. A non-zero
813 : : //! vector result indicates a leak, e.g., a hole in the partition (covered by
814 : : //! the faces of the face adjacency communication map), which indicates an
815 : : //! error upstream in the code that sets up the face communication data
816 : : //! structures.
817 : : //! \note Compared to tk::leakyPartition() this function performs the leak-test
818 : : //! on the face geometry data structure enlarged by ghost faces on this
819 : : //! partition by computing a discrete surface integral considering the
820 : : //! physical and chare boundary faces, which should be equal to zero for a
821 : : //! closed domain.
822 : : //! \return True if our chare face adjacency leaks.
823 : : // *****************************************************************************
824 : : {
825 : : // Storage for surface integral over our chunk of the adjacency
826 : 1099 : std::array< tk::real, 3 > s{{ 0.0, 0.0, 0.0 }};
827 : :
828 : : // physical boundary faces
829 [ + - ][ + + ]: 123789 : for (std::size_t f=0; f<m_fd.Nbfac(); ++f) {
830 [ + - ][ + - ]: 122690 : s[0] += m_geoFace(f,0) * m_geoFace(f,1);
831 [ + - ][ + - ]: 122690 : s[1] += m_geoFace(f,0) * m_geoFace(f,2);
832 [ + - ][ + - ]: 122690 : s[2] += m_geoFace(f,0) * m_geoFace(f,3);
833 : : }
834 : :
835 : : // chare-boundary faces
836 [ + + ]: 38755 : for (std::size_t f=m_fd.Nipfac(); f<m_fd.Esuf().size()/2; ++f) {
837 [ + - ][ + - ]: 37656 : s[0] += m_geoFace(f,0) * m_geoFace(f,1);
838 [ + - ][ + - ]: 37656 : s[1] += m_geoFace(f,0) * m_geoFace(f,2);
839 [ + - ][ + - ]: 37656 : s[2] += m_geoFace(f,0) * m_geoFace(f,3);
840 : : }
841 : :
842 : 1099 : auto eps = std::numeric_limits< tk::real >::epsilon() * 100;
843 [ + - ][ + - ]: 1099 : return std::abs(s[0]) > eps || std::abs(s[1]) > eps || std::abs(s[2]) > eps;
[ - + ]
844 : : }
845 : :
846 : : bool
847 : 1099 : Ghosts::faceMatch()
848 : : // *****************************************************************************
849 : : // Check if esuf of chare-boundary faces matches
850 : : //! \details This function checks each chare-boundary esuf entry for the left
851 : : //! and right elements. Then, it tries to match all vertices of these
852 : : //! elements. Exactly three of these vertices must match if the esuf entry
853 : : //! has been updated correctly at chare-boundaries.
854 : : //! \return True if chare-boundary faces match.
855 : : // *****************************************************************************
856 : : {
857 : 1099 : const auto& esuf = m_fd.Esuf();
858 : 1099 : bool match(true);
859 : :
860 : 1099 : auto eps = std::numeric_limits< tk::real >::epsilon() * 100;
861 : :
862 [ + + ]: 38755 : for (auto f=m_fd.Nipfac(); f<esuf.size()/2; ++f)
863 : : {
864 : 37656 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
865 : 37656 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
866 : :
867 : 37656 : std::size_t count = 0;
868 : :
869 [ + + ]: 188280 : for (std::size_t i=0; i<4; ++i)
870 : : {
871 : 150624 : auto ip = m_inpoel[4*el+i];
872 [ + + ]: 753120 : for (std::size_t j=0; j<4; ++j)
873 : : {
874 : 602496 : auto jp = m_inpoel[4*er+j];
875 : 602496 : auto xdiff = std::abs( m_coord[0][ip] - m_coord[0][jp] );
876 : 602496 : auto ydiff = std::abs( m_coord[1][ip] - m_coord[1][jp] );
877 : 602496 : auto zdiff = std::abs( m_coord[2][ip] - m_coord[2][jp] );
878 : :
879 [ + + ][ + + ]: 602496 : if ( xdiff<=eps && ydiff<=eps && zdiff<=eps ) ++count;
[ + + ]
880 : : }
881 : : }
882 : :
883 [ + - ][ + - ]: 37656 : match = (match && count == 3);
884 : : }
885 : :
886 : 1099 : return match;
887 : : }
888 : :
889 : : bool
890 : 1061 : Ghosts::receivedChBndFaces()
891 : : // *****************************************************************************
892 : : // Verify that all chare-boundary faces have been received
893 : : //! \return True if all chare-boundary faces have been received
894 : : // *****************************************************************************
895 : : {
896 [ + - ]: 1061 : const auto& lid = Disc()->Lid();
897 : 2122 : tk::UnsMesh::FaceSet recvBndFace;
898 : :
899 : : // Collect chare-boundary faces that have been received and expected
900 [ + + ]: 6387 : for (const auto& c : m_bndFace)
901 [ + + ]: 42982 : for (const auto& f : c.second)
902 [ + - ][ + - ]: 37656 : if (m_expChBndFace.find(f.first) != end(m_expChBndFace))
903 [ + - ]: 37656 : recvBndFace.insert(f.first);
904 : :
905 : : // Collect info on expected but not received faces
906 [ + - ]: 2122 : std::stringstream msg;
907 [ + + ]: 38717 : for (const auto& f : m_expChBndFace)
908 [ + - ][ - + ]: 37656 : if (recvBndFace.find(f) == end(recvBndFace)) {
909 : 0 : const auto& x = m_coord[0];
910 : 0 : const auto& y = m_coord[1];
911 : 0 : const auto& z = m_coord[2];
912 [ - - ]: 0 : auto A = tk::cref_find( lid, f[0] );
913 [ - - ]: 0 : auto B = tk::cref_find( lid, f[1] );
914 [ - - ]: 0 : auto C = tk::cref_find( lid, f[2] );
915 [ - - ][ - - ]: 0 : msg << '{' << A << ',' << B << ',' << C << "}:("
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
916 [ - - ][ - - ]: 0 : << x[A] << ',' << y[A] << ',' << z[A] << ' '
[ - - ][ - - ]
[ - - ][ - - ]
917 [ - - ][ - - ]: 0 : << x[B] << ',' << y[B] << ',' << z[B] << ' '
[ - - ][ - - ]
[ - - ][ - - ]
918 [ - - ][ - - ]: 0 : << x[C] << ',' << y[C] << ',' << z[C] << ") ";
[ - - ][ - - ]
[ - - ][ - - ]
919 : : }
920 : :
921 : 1061 : tk::destroy( m_expChBndFace );
922 : :
923 : : // Error out with info on missing faces
924 [ + - ]: 1061 : auto s = msg.str();
925 [ - + ]: 1061 : if (!s.empty()) {
926 [ - - ][ - - ]: 0 : Throw( "Ghosts chare " + std::to_string(thisIndex) +
[ - - ][ - - ]
[ - - ][ - - ]
927 : : " missing face(s) {local node ids} (node coords): " + s );
928 : : } else {
929 : 2122 : return true;
930 : : }
931 : : }
932 : :
933 : : int
934 : 215680 : Ghosts::findchare( const tk::UnsMesh::Face& t )
935 : : // *****************************************************************************
936 : : // Find any chare for face (given by 3 global node IDs)
937 : : //! \param[in] t Face given by three global node IDs
938 : : //! \return Chare ID if found among any of the chares we communicate along
939 : : //! faces with, -1 if the face cannot be found.
940 : : // *****************************************************************************
941 : : {
942 [ + + ]: 800284 : for (const auto& cf : m_bndFace)
943 : : // cppcheck-suppress useStlAlgorithm
944 [ + - ][ + + ]: 659916 : if (cf.second.find(t) != end(cf.second))
945 : 75312 : return cf.first;
946 : 140368 : return -1;
947 : : }
948 : :
949 : : std::size_t
950 : 37656 : Ghosts::nodetripletMatch(
951 : : const std::array< std::size_t, 2 >& id,
952 : : const tk::UnsMesh::Face& t )
953 : : // *****************************************************************************
954 : : // Check if entries in inpoel, inpofa and node-triplet are consistent
955 : : //! \param[in] id Local face and (inner) tet id adjacent to it
956 : : //! \param[in] t node-triplet associated with the chare boundary face
957 : : //! \return number of nodes in inpoel that matched with t and inpofa
958 : : // *****************************************************************************
959 : : {
960 : 37656 : const auto& esuf = m_fd.Esuf();
961 : 37656 : const auto& inpofa = m_fd.Inpofa();
962 : 37656 : const auto& lid = Disc()->Lid();
963 : :
964 : 37656 : std::size_t counter = 0;
965 [ + + ]: 188280 : for (std::size_t k=0; k<4; ++k)
966 : : {
967 : 150624 : auto el = esuf[ 2*id[0] ];
968 : 150624 : auto ip = m_inpoel[ 4*static_cast< std::size_t >( el )+k ];
969 [ - + ][ - - ]: 150624 : Assert( el == static_cast< int >( id[1] ), "Mismatch in id and esuf" );
[ - - ][ - - ]
970 [ + + ]: 602496 : for (std::size_t j=0; j<3; ++j)
971 : : {
972 : 451872 : auto jp = tk::cref_find( lid, t[j] );
973 : 451872 : auto fp = inpofa[ 3*id[0]+(2-j) ];
974 [ + + ][ + - ]: 451872 : if (ip == jp && ip == fp) ++counter;
975 : : }
976 : : }
977 : :
978 : 37656 : return counter;
979 : : }
980 : :
981 : : void
982 : 37656 : Ghosts::addEsuf(
983 : : const std::array< std::size_t, 2 >& id,
984 : : std::size_t ghostid )
985 : : // *****************************************************************************
986 : : // Fill elements surrounding a face along chare boundary
987 : : //! \param[in] id Local face and (inner) tet id adjacent to it
988 : : //! \param[in] ghostid Local ID for ghost tet
989 : : //! \details This function extends and fills in the elements surrounding faces
990 : : //! data structure (esuf) so that the left and right element id is filled
991 : : //! in correctly on chare boundaries to contain the correct inner tet id and
992 : : //! the local tet id for the outer (ghost) tet, both adjacent to the given
993 : : //! chare-face boundary. Prior to this function, this data structure does not
994 : : //! have yet face-element connectivity adjacent to chare-boundary faces, only
995 : : //! for physical boundaries and internal faces that are not on the chare
996 : : //! boundary (this latter purely as a result of mesh partitioning). The remote
997 : : //! element id of the ghost is stored in a location that is local to our own
998 : : //! esuf. The face numbering is such that esuf stores the element-face
999 : : //! connectivity first for the physical-boundary faces, followed by that of
1000 : : //! the internal faces, followed by the chare-boundary faces. As a result,
1001 : : //! esuf can be used by physics algorithms in exactly the same way as would be
1002 : : //! used in serial. In serial, of course, this data structure is not extended
1003 : : //! at the end by the chare-boundaries.
1004 : : // *****************************************************************************
1005 : : {
1006 : 37656 : auto& esuf = m_fd.Esuf();
1007 : :
1008 [ - + ][ - - ]: 37656 : Assert( 2*id[0]+1 < esuf.size(), "Indexing out of esuf" );
[ - - ][ - - ]
1009 : :
1010 : : // put in inner tet id
1011 [ + - ][ + - ]: 37656 : Assert( esuf[ 2*id[0] ] == -2 && esuf[ 2*id[0]+1 ] == -2, "Updating esuf at "
[ - - ][ - - ]
[ - - ]
1012 : : "wrong location instead of chare-boundary" );
1013 : 37656 : esuf[ 2*id[0]+0 ] = static_cast< int >( id[1] );
1014 : : // put in local id for outer/ghost tet
1015 : 37656 : esuf[ 2*id[0]+1 ] = static_cast< int >( ghostid );
1016 : 37656 : }
1017 : :
1018 : : void
1019 : 37656 : Ghosts::addEsuel(
1020 : : const std::array< std::size_t, 2 >& id,
1021 : : std::size_t ghostid,
1022 : : const tk::UnsMesh::Face& t )
1023 : : // *****************************************************************************
1024 : : // Fill elements surrounding a element along chare boundary
1025 : : //! \param[in] id Local face and (inner) tet id adjacent to it
1026 : : //! \param[in] ghostid Local ID for ghost tet
1027 : : //! \param[in] t node-triplet associated with the chare boundary face
1028 : : //! \details This function updates the elements surrounding element (esuel) data
1029 : : // structure for the (inner) tets adjacent to the chare-boundaries. It fills
1030 : : // esuel of this inner tet with the local tet-id that has been assigned to
1031 : : // the outer ghost tet in Ghosts::comGhost in place of the -1 before.
1032 : : // *****************************************************************************
1033 : : {
1034 [ + - ]: 37656 : const auto& lid = Disc()->Lid();
1035 : 37656 : [[maybe_unused]] const auto& esuf = m_fd.Esuf();
1036 : 37656 : auto& esuel = m_fd.Esuel();
1037 : :
1038 : : std::array< tk::UnsMesh::Face, 4 > face;
1039 [ + + ]: 188280 : for (std::size_t f = 0; f<4; ++f)
1040 [ + + ]: 602496 : for (std::size_t i = 0; i<3; ++i)
1041 : 451872 : face[f][i] = m_inpoel[ id[1]*4 + tk::lpofa[f][i] ];
1042 : :
1043 [ + - ]: 37656 : tk::UnsMesh::Face tl{{ tk::cref_find( lid, t[0] ),
1044 [ + - ]: 37656 : tk::cref_find( lid, t[1] ),
1045 [ + - ]: 75312 : tk::cref_find( lid, t[2] ) }};
1046 : :
1047 : 37656 : std::size_t i(0), nmatch(0);
1048 [ + + ]: 188280 : for (const auto& f : face) {
1049 [ + - ][ + + ]: 150624 : if (tk::UnsMesh::Eq< 3 >()( tl, f )) {
1050 [ - + ][ - - ]: 37656 : Assert( esuel[ id[1]*4 + i ] == -1, "Incorrect boundary element found in "
[ - - ][ - - ]
1051 : : "esuel");
1052 : 37656 : esuel[ id[1]*4 + i ] = static_cast<int>(ghostid);
1053 : 37656 : ++nmatch;
1054 [ - + ][ - - ]: 37656 : Assert( esuel[ id[1]*4 + i ] == esuf[ 2*id[0]+1 ], "Incorrect boundary "
[ - - ][ - - ]
1055 : : "element entered in esuel" );
1056 [ - + ][ - - ]: 37656 : Assert( static_cast<int>(id[1]) == esuf[ 2*id[0]+0 ], "Boundary "
[ - - ][ - - ]
1057 : : "element entered in incorrect esuel location" );
1058 : : }
1059 : 150624 : ++i;
1060 : : }
1061 : :
1062 : : // ensure that exactly one face matched
1063 [ - + ][ - - ]: 37656 : Assert( nmatch == 1, "Incorrect number of node-triplets (faces) matched for "
[ - - ][ - - ]
[ - - ]
1064 : : "updating esuel; matching faces = "+ std::to_string(nmatch) );
1065 : 37656 : }
1066 : :
1067 : : void
1068 : 37656 : Ghosts::addGeoFace(
1069 : : const tk::UnsMesh::Face& t,
1070 : : const std::array< std::size_t, 2 >& id )
1071 : : // *****************************************************************************
1072 : : // Fill face-geometry data along chare boundary
1073 : : //! \param[in] t Face (given by 3 global node IDs) on the chare boundary
1074 : : //! \param[in] id Local face and (inner) tet id adjacent to face t
1075 : : //! \details This function fills in the face geometry data along a chare
1076 : : //! boundary.
1077 : : // *****************************************************************************
1078 : : {
1079 [ + - ]: 37656 : const auto& lid = Disc()->Lid();
1080 : :
1081 : : // get global node IDs reversing order to get outward-pointing normal
1082 [ + - ]: 37656 : auto A = tk::cref_find( lid, t[2] );
1083 [ + - ]: 37656 : auto B = tk::cref_find( lid, t[1] );
1084 [ + - ]: 37656 : auto C = tk::cref_find( lid, t[0] );
1085 : 112968 : auto geochf = tk::geoFaceTri( {{m_coord[0][A], m_coord[0][B], m_coord[0][C]}},
1086 : 112968 : {{m_coord[1][A], m_coord[1][B], m_coord[1][C]}},
1087 [ + - ]: 301248 : {{m_coord[2][A], m_coord[2][B], m_coord[2][C]}} );
1088 : :
1089 [ + + ]: 301248 : for (std::size_t i=0; i<7; ++i)
1090 [ + - ][ + - ]: 263592 : m_geoFace(id[0],i) = geochf(0,i);
1091 : 37656 : }
1092 : :
1093 : : #include "NoWarning/ghosts.def.h"
|
