Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/Partitioner.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 Charm++ chare partitioner nodegroup used to perform mesh
9 : : partitioning
10 : : \details Charm++ chare partitioner nodegroup used to perform mesh read and
11 : : partitioning, one worker per compute node.
12 : : */
13 : : // *****************************************************************************
14 : :
15 : : #include <numeric>
16 : :
17 : : #include "NoWarning/kokkos_cr.hpp"
18 : :
19 : : #include "Partitioner.hpp"
20 : : #include "DerivedData.hpp"
21 : : #include "Reorder.hpp"
22 : : #include "MeshReader.hpp"
23 : : #include "CGPDE.hpp"
24 : : #include "DGPDE.hpp"
25 : : #include "Inciter/Options/Scheme.hpp"
26 : : #include "UnsMesh.hpp"
27 : : #include "Vector.hpp"
28 : : #include "ContainerUtil.hpp"
29 : : #include "Callback.hpp"
30 : :
31 : : namespace inciter {
32 : :
33 : : extern ctr::InputDeck g_inputdeck;
34 : :
35 : : } // inciter::
36 : :
37 : : using inciter::Partitioner;
38 : :
39 : 527 : Partitioner::Partitioner(
40 : : std::size_t meshid,
41 : : const std::string& filename,
42 : : const tk::PartitionerCallback& cbp,
43 : : const tk::RefinerCallback& cbr,
44 : : const tk::SorterCallback& cbs,
45 : : const CProxy_Transporter& host,
46 : : const CProxy_Refiner& refiner,
47 : : const CProxy_Sorter& sorter,
48 : : const tk::CProxy_MeshWriter& meshwriter,
49 : : const std::vector< Scheme >& scheme,
50 : : const std::map< int, std::vector< std::size_t > >& bface,
51 : : const std::map< int, std::vector< std::size_t > >& faces,
52 : 527 : const std::map< int, std::vector< std::size_t > >& bnode ) :
53 : : m_meshid( meshid ),
54 : : m_cbp( cbp ),
55 : : m_cbr( cbr ),
56 : : m_cbs( cbs ),
57 : : m_host( host ),
58 : : m_refiner( refiner ),
59 : : m_sorter( sorter ),
60 : : m_meshwriter( meshwriter ),
61 : : m_scheme( scheme ),
62 : : m_ginpoel(),
63 : : m_coord(),
64 : : m_inpoel(),
65 : : m_lid(),
66 : : m_elemBlockId(),
67 : : m_ndist( 0 ),
68 : : m_nchare( 0 ),
69 : : m_nface(),
70 : : m_chinpoel(),
71 : : m_chcoordmap(),
72 : : m_chbface(),
73 : : m_chtriinpoel(),
74 : : m_chbnode(),
75 : : m_chelemblockid(),
76 : : m_bface( bface ),
77 [ + - ][ + - ]: 1581 : m_bnode( bnode )
[ + - ][ + - ]
78 : : // *****************************************************************************
79 : : // Constructor
80 : : //! \param[in] meshid Mesh ID
81 : : //! \param[in] filename Input mesh filename to read from
82 : : //! \param[in] cbp Charm++ callbacks for Partitioner
83 : : //! \param[in] cbr Charm++ callbacks for Refiner
84 : : //! \param[in] cbs Charm++ callbacks for Sorter
85 : : //! \param[in] host Host Charm++ proxy we are being called from
86 : : //! \param[in] refiner Mesh refiner proxy
87 : : //! \param[in] sorter Mesh reordering (sorter) proxy
88 : : //! \param[in] meshwriter Mesh writer proxy
89 : : //! \param[in] scheme Discretization scheme
90 : : //! \param[in] bface File-internal elem ids of side sets (whole mesh)
91 : : //! \param[in] faces Elem-relative face ids of side sets (whole mesh)
92 : : //! \param[in] bnode Node lists of side sets (whole mesh)
93 : : // *****************************************************************************
94 : : {
95 : : // The following call has to be made on all MPI ranks immediately after
96 : : // initializing MPI. This has to be done from a Charm++ nodegroup, since there
97 : : // exists only one rank per node on a nodegroup, in SMP mode.
98 : : // see https://github.com/trilinos/Trilinos/issues/11197#issuecomment-1301325163
99 [ + + ][ + - ]: 527 : if (m_meshid == 0) Kokkos::initialize();
100 : :
101 : : // Create mesh reader
102 [ + - ]: 527 : tk::MeshReader mr( filename );
103 : :
104 : : // Read this compute node's chunk of the mesh (graph and coords) from file
105 : : std::vector< std::size_t > triinpoel;
106 [ + - ]: 527 : mr.readMeshPart( m_ginpoel, m_inpoel, triinpoel, m_lid, m_coord,
107 : : m_elemBlockId, CkNumNodes(), CkMyNode() );
108 : :
109 : : // Check if mesh rotation specified
110 : : auto mesh_orientation =
111 [ + - ]: 527 : g_inputdeck.get< tag::mesh >()[meshid].get< tag::orientation >();
112 : : bool rotate_mesh(false);
113 [ + + ]: 2108 : for (std::size_t i=0; i<3; ++i) {
114 [ + - ]: 1581 : if (std::abs(mesh_orientation[i]) > 1e-8) {
115 : : rotate_mesh = true;
116 : : break;
117 : : }
118 : : }
119 : :
120 : : // Rotate mesh if specified
121 [ - + ]: 527 : if (rotate_mesh) {
122 : : auto mesh_cm =
123 [ - - ]: 0 : g_inputdeck.get< tag::mesh >()[meshid].get< tag::center_of_mass >();
124 : : auto& x = m_coord[0];
125 : : auto& y = m_coord[1];
126 : : auto& z = m_coord[2];
127 [ - - ]: 0 : for (std::size_t i=0; i<x.size(); ++i) {
128 : : std::array< tk::real, 3 >
129 : 0 : point{{ x[i]-mesh_cm[0], y[i]-mesh_cm[1], z[i]-mesh_cm[2] }};
130 : 0 : tk::rotatePoint(
131 : 0 : {{ mesh_orientation[0], mesh_orientation[1], mesh_orientation[2] }},
132 : : point );
133 : 0 : x[i] = point[0]+mesh_cm[0];
134 : 0 : y[i] = point[1]+mesh_cm[1];
135 : 0 : z[i] = point[2]+mesh_cm[2];
136 : : }
137 : : }
138 : :
139 : : // Check if mesh location specified
140 : : auto mesh_location =
141 [ + - ]: 527 : g_inputdeck.get< tag::mesh >()[meshid].get< tag::location >();
142 : : bool relocate_mesh(false);
143 [ + + ]: 2108 : for (std::size_t i=0; i<3; ++i) {
144 [ + - ]: 1581 : if (std::abs(mesh_location[i]) > 1e-8) {
145 : : relocate_mesh = true;
146 : : break;
147 : : }
148 : : }
149 : :
150 : : // Relocate mesh if specified
151 [ - + ]: 527 : if (relocate_mesh) {
152 : : auto& x = m_coord[0];
153 : : auto& y = m_coord[1];
154 : : auto& z = m_coord[2];
155 [ - - ]: 0 : for (std::size_t i=0; i<x.size(); ++i) {
156 : 0 : std::array< tk::real, 3 > point{{ x[i], y[i], z[i] }};
157 : : // negative values due to how tk::movePoint() is interpreted
158 : 0 : tk::movePoint(
159 : 0 : {{ -mesh_location[0], -mesh_location[1], -mesh_location[2] }}, point );
160 : 0 : x[i] = point[0];
161 : 0 : y[i] = point[1];
162 : 0 : z[i] = point[2];
163 : : }
164 : : }
165 : :
166 : : // Compute triangle connectivity for side sets, reduce boundary face for side
167 : : // sets to this compute node only and to compute-node-local face ids
168 : 527 : m_triinpoel = mr.triinpoel( m_bface, faces, m_ginpoel, triinpoel );
169 : :
170 : : // Reduce boundary node lists (global ids) for side sets to this compute node
171 : : // only
172 [ + - ]: 527 : ownBndNodes( m_lid, m_bnode );
173 : :
174 : : // Sum number of cells across distributed mesh
175 [ + - ][ - - ]: 527 : std::vector< std::size_t > meshdata{ meshid, m_ginpoel.size()/4 };
176 [ + - ]: 527 : contribute( meshdata, CkReduction::sum_ulong, m_cbp.get< tag::load >() );
177 : 527 : }
178 : :
179 : : void
180 : 527 : Partitioner::ownBndNodes(
181 : : const std::unordered_map< std::size_t, std::size_t >& lid,
182 : : std::map< int, std::vector< std::size_t > >& bnode )
183 : : // *****************************************************************************
184 : : // Keep only those nodes for side sets that reside on this compute node
185 : : //! \param[in] lid Global->local node IDs of elements of this compute node's
186 : : //! mesh chunk
187 : : //! \param[in,out] bnode Global ids of nodes for side sets for whole mesh
188 : : //! \details This function overwrites the input boundary node lists map with the
189 : : //! nodes that reside on the caller compute node.
190 : : // *****************************************************************************
191 : : {
192 : : std::map< int, std::vector< std::size_t > > bnode_own;
193 : :
194 [ + + ]: 1755 : for (const auto& [ setid, nodes ] : bnode) {
195 [ + - ]: 1228 : auto& b = bnode_own[ setid ];
196 [ + + ]: 112027 : for (auto n : nodes) {
197 : : auto i = lid.find( n );
198 [ + + ][ + - ]: 110799 : if (i != end(lid)) b.push_back( n );
199 : : }
200 [ + + ]: 1228 : if (b.empty()) bnode_own.erase( setid );
201 : : }
202 : :
203 : : bnode = std::move(bnode_own);
204 : 527 : }
205 : :
206 : : void
207 : 527 : Partitioner::partition( int nchare )
208 : : // *****************************************************************************
209 : : // Partition the computational mesh into a number of chares
210 : : //! \param[in] nchare Number of parts the mesh will be partitioned into
211 : : //! \details This function calls the mesh partitioner to partition the mesh. The
212 : : //! number of partitions equals the number nchare argument which must be no
213 : : //! lower than the number of compute nodes.
214 : : // *****************************************************************************
215 : : {
216 : : Assert( nchare >= CkNumNodes(), "Number of chares must not be lower than the "
217 : : "number of compute nodes" );
218 : :
219 : : // Generate element IDs for Zoltan
220 : 527 : std::vector< long long > gelemid( m_ginpoel.size()/4 );
221 : : std::iota( begin(gelemid), end(gelemid), 0 );
222 : :
223 : 527 : m_nchare = nchare;
224 : 527 : const auto alg = g_inputdeck.get< tag::partitioning >();
225 : : const auto che = tk::zoltan::geomPartMesh( alg,
226 [ + - ][ - - ]: 527 : centroids( m_inpoel, m_coord ),
227 : : gelemid,
228 [ + - ]: 527 : nchare );
229 : :
230 [ - + ][ - - ]: 527 : if ( g_inputdeck.get< tag::cmd, tag::feedback >() ) m_host.pepartitioned();
231 : :
232 [ + - ]: 527 : contribute( sizeof(std::size_t), &m_meshid, CkReduction::nop,
233 : : m_cbp.get< tag::partitioned >() );
234 : :
235 : : Assert( che.size() == gelemid.size(), "Size of ownership array (chare ID "
236 : : "of elements) after mesh partitioning does not equal the number of "
237 : : "mesh graph elements" );
238 : :
239 : : // Categorize mesh elements (given by their gobal node IDs) by target chare
240 : : // and distribute to their compute nodes based on mesh partitioning.
241 [ + - ][ + - ]: 1054 : distribute( categorize( che ) );
[ + - ][ - - ]
242 : 527 : }
243 : :
244 : : void
245 : 1992 : Partitioner::addMesh(
246 : : int fromnode,
247 : : const std::unordered_map< int, // chare id
248 : : std::tuple<
249 : : std::vector< std::size_t >, // tet connectivity
250 : : tk::UnsMesh::CoordMap, // node coords
251 : : std::unordered_map< int, std::vector< std::size_t > >, // bface conn
252 : : std::unordered_map< int, std::vector< std::size_t > >, // bnodes
253 : : std::vector< std::size_t > // elem-blocks
254 : : > >& chmesh )
255 : : // *****************************************************************************
256 : : // Receive mesh associated to chares we own after refinement
257 : : //! \param[in] fromnode Compute node call coming from
258 : : //! \param[in] chmesh Map associating mesh connectivities to global node ids
259 : : //! and node coordinates for mesh chunks we are assigned by the partitioner
260 : : // *****************************************************************************
261 : : {
262 : : // Store mesh connectivity and global node coordinates categorized by chares.
263 : : // The send side also writes to the data written here, so concat.
264 [ + + ]: 7642 : for (const auto& [ chareid, chunk ] : chmesh) {
265 : : Assert( node(chareid) == CkMyNode(), "Compute node "
266 : : + std::to_string(CkMyNode()) +
267 : : " received a mesh whose chare it does not own" );
268 : : // Store domain element (tetrahedron) connectivity
269 : : const auto& inpoel = std::get< 0 >( chunk );
270 : : auto& inp = m_chinpoel[ chareid ]; // will store tetrahedron connectivity
271 : 5650 : inp.insert( end(inp), begin(inpoel), end(inpoel) );
272 : : // Store own mesh block id and associated tet ids
273 : : const auto& elblock = std::get<4>( chunk );
274 : : auto& cheb = m_chelemblockid[ chareid ];
275 : 5650 : cheb.insert( end(cheb), begin(elblock), end(elblock) );
276 : : // Store mesh node coordinates associated to global node IDs
277 : : const auto& coord = std::get< 1 >( chunk );
278 : : Assert( tk::uniquecopy(inpoel).size() == coord.size(), "Size mismatch" );
279 : : auto& chcm = m_chcoordmap[ chareid ]; // will store node coordinates
280 : : chcm.insert( begin(coord), end(coord) ); // concatenate node coords
281 : : // Store boundary side set id + face ids + face connectivities
282 : : const auto& bconn = std::get< 2 >( chunk );
283 : : auto& bface = m_chbface[ chareid ]; // for side set id + boundary face ids
284 : 5650 : auto& t = m_chtriinpoel[ chareid ]; // for boundary face connectivity
285 : : auto& f = m_nface[ chareid ]; // use counter for chare
286 [ + + ]: 13597 : for (const auto& [ setid, faceids ] : bconn) {
287 : 7947 : auto& b = bface[ setid ];
288 [ + + ]: 64889 : for (std::size_t i=0; i<faceids.size()/3; ++i) {
289 : 56942 : b.push_back( f++ );
290 : 56942 : t.push_back( faceids[i*3+0] );
291 : 56942 : t.push_back( faceids[i*3+1] );
292 : 56942 : t.push_back( faceids[i*3+2] );
293 : : }
294 : : }
295 : : // Store boundary side set id + node lists
296 : : const auto& bnode = std::get< 3 >( chunk );
297 : : auto& nodes = m_chbnode[ chareid ]; // for side set id + boundary nodes
298 [ + + ]: 10660 : for (const auto& [ setid, bnodes ] : bnode) {
299 : 5010 : auto& b = nodes[ setid ];
300 : 5010 : b.insert( end(b), begin(bnodes), end(bnodes) );
301 : : }
302 : : }
303 : :
304 [ + - ]: 1992 : thisProxy[ fromnode ].recvMesh();
305 : 1992 : }
306 : :
307 : : int
308 : 5650 : Partitioner::node( int id ) const
309 : : // *****************************************************************************
310 : : // Return nodegroup id for chare id
311 : : //! \param[in] id Chare id
312 : : //! \return Nodegroup that creates the chare
313 : : //! \details This is computed based on a simple contiguous linear
314 : : //! distribution of chare ids to compute nodes.
315 : : // *****************************************************************************
316 : : {
317 : : Assert( m_nchare > 0, "Number of chares must be a positive number" );
318 : 5650 : auto p = id / (m_nchare / CkNumNodes());
319 [ + + ]: 5650 : if (p >= CkNumNodes()) p = CkNumNodes()-1;
320 : : Assert( p < CkNumNodes(), "Assigning to nonexistent node" );
321 : 5650 : return p;
322 : : }
323 : :
324 : : void
325 : 1992 : Partitioner::recvMesh()
326 : : // *****************************************************************************
327 : : // Acknowledge received mesh chunk and its nodes after mesh refinement
328 : : // *****************************************************************************
329 : : {
330 [ + + ]: 1992 : if (--m_ndist == 0) {
331 [ - + ]: 450 : if (g_inputdeck.get< tag::cmd, tag::feedback >()) m_host.pedistributed();
332 : 450 : contribute( sizeof(std::size_t), &m_meshid, CkReduction::nop,
333 : : m_cbp.get< tag::distributed >() );
334 : : }
335 : 1992 : }
336 : :
337 : : void
338 : 527 : Partitioner::refine()
339 : : // *****************************************************************************
340 : : // Optionally start refining the mesh
341 : : // *****************************************************************************
342 : : {
343 [ + - ]: 527 : auto dist = distribution( m_nchare );
344 : :
345 : : std::size_t error = 0;
346 [ + - ]: 527 : if (m_chinpoel.size() < static_cast<std::size_t>(dist[1])) {
347 : :
348 : : error = 1;
349 : :
350 : : } else {
351 : :
352 [ + + ]: 2235 : for (int c=0; c<dist[1]; ++c) {
353 : : // compute chare ID
354 : 1708 : auto cid = CkMyNode() * dist[0] + c;
355 : :
356 : : Assert(tk::cref_find(m_chinpoel,cid).size()/4 ==
357 : : tk::cref_find(m_chelemblockid,cid).size(),
358 : : "incorrect number of elements in elem-block id vector");
359 : : // create refiner Charm++ chare array element using dynamic insertion
360 : 3416 : m_refiner[ cid ].insert( m_meshid,
361 : 1708 : m_host,
362 : 1708 : m_sorter,
363 : 1708 : m_meshwriter,
364 : 1708 : m_scheme,
365 : 1708 : m_cbr,
366 [ + - ]: 1708 : m_cbs,
367 : : tk::cref_find(m_chinpoel,cid),
368 : : tk::cref_find(m_chcoordmap,cid),
369 : : tk::cref_find(m_chbface,cid),
370 : : tk::cref_find(m_chtriinpoel,cid),
371 : : tk::cref_find(m_chbnode,cid),
372 : : tk::cref_find(m_chelemblockid,cid),
373 : : m_nchare );
374 : : }
375 : :
376 : : }
377 : :
378 : : tk::destroy( m_ginpoel );
379 : 527 : tk::destroy( m_coord );
380 : : tk::destroy( m_inpoel );
381 : 527 : tk::destroy( m_lid );
382 : 527 : tk::destroy( m_elemBlockId );
383 : 527 : tk::destroy( m_nface );
384 : 527 : tk::destroy( m_nodech );
385 : 527 : tk::destroy( m_linnodes );
386 : 527 : tk::destroy( m_chinpoel );
387 : 527 : tk::destroy( m_chcoordmap );
388 : 527 : tk::destroy( m_chbface );
389 : 527 : tk::destroy( m_chtriinpoel );
390 : 527 : tk::destroy( m_chbnode );
391 : 527 : tk::destroy( m_chelemblockid );
392 : 527 : tk::destroy( m_bnodechares );
393 : 527 : tk::destroy( m_bface );
394 : : tk::destroy( m_triinpoel );
395 : 527 : tk::destroy( m_bnode );
396 : :
397 : 527 : std::vector< std::size_t > meshdata{ m_meshid, error };
398 [ + - ]: 527 : contribute( meshdata, CkReduction::max_ulong, m_cbp.get<tag::refinserted>() );
399 : 527 : }
400 : :
401 : : std::array< std::vector< tk::real >, 3 >
402 : 527 : Partitioner::centroids( const std::vector< std::size_t >& inpoel,
403 : : const tk::UnsMesh::Coords& coord )
404 : : // *****************************************************************************
405 : : // Compute element centroid coordinates
406 : : //! \param[in] inpoel Mesh connectivity with local ids
407 : : //! \param[in] coord Node coordinates
408 : : //! \return Centroids for all cells on this compute node
409 : : // *****************************************************************************
410 : : {
411 : : Assert( tk::uniquecopy(inpoel).size() == coord[0].size(), "Size mismatch" );
412 : :
413 : : const auto& x = coord[0];
414 : : const auto& y = coord[1];
415 : : const auto& z = coord[2];
416 : :
417 : : // Make room for element centroid coordinates
418 : : std::array< std::vector< tk::real >, 3 > cent;
419 : : auto& cx = cent[0];
420 : : auto& cy = cent[1];
421 : : auto& cz = cent[2];
422 [ + - ]: 527 : auto num = inpoel.size()/4;
423 [ + - ]: 527 : cx.resize( num );
424 [ + - ]: 527 : cy.resize( num );
425 [ + - ]: 527 : cz.resize( num );
426 : :
427 : : // Compute element centroids for mesh passed in
428 [ + + ]: 310333 : for (std::size_t e=0; e<num; ++e) {
429 : 309806 : auto A = inpoel[e*4+0];
430 : 309806 : auto B = inpoel[e*4+1];
431 : 309806 : auto C = inpoel[e*4+2];
432 : 309806 : auto D = inpoel[e*4+3];
433 : 309806 : cx[e] = (x[A] + x[B] + x[C] + x[D]) / 4.0;
434 : 309806 : cy[e] = (y[A] + y[B] + y[C] + y[D]) / 4.0;
435 : 309806 : cz[e] = (z[A] + z[B] + z[C] + z[D]) / 4.0;
436 : : }
437 : :
438 : 527 : return cent;
439 : : }
440 : :
441 : : std::unordered_map< int, Partitioner::MeshData >
442 : 527 : Partitioner::categorize( const std::vector< std::size_t >& target ) const
443 : : // *****************************************************************************
444 : : // Categorize mesh data by target
445 : : //! \param[in] target Target chares of mesh elements, size: number of
446 : : //! elements in the chunk of the mesh graph on this compute node.
447 : : //! \return Vector of global mesh node ids connecting elements owned by each
448 : : //! target chare.
449 : : // *****************************************************************************
450 : : {
451 : : Assert( target.size() == m_ginpoel.size()/4, "Size mismatch");
452 : :
453 : : using Face = tk::UnsMesh::Face;
454 : :
455 : : // Build hash map associating side set id to boundary faces
456 : : std::unordered_map< Face, int,
457 : : tk::UnsMesh::Hash<3>, tk::UnsMesh::Eq<3> > faceside;
458 [ + + ]: 3245 : for (const auto& [ setid, faceids ] : m_bface)
459 [ + + ]: 137200 : for (auto f : faceids)
460 : 134482 : faceside[ {{ m_triinpoel[f*3+0],
461 [ + - ]: 134482 : m_triinpoel[f*3+1],
462 [ + - ]: 134482 : m_triinpoel[f*3+2] }} ] = setid;
463 : :
464 : : // Build hash map associating side set ids to boundary nodes
465 : : std::unordered_map< std::size_t, std::unordered_set< int > > nodeside;
466 [ + + ]: 1753 : for (const auto& [ setid, nodes ] : m_bnode)
467 [ + + ][ + - ]: 79644 : for (auto n : nodes)
468 : : nodeside[ n ].insert( setid );
469 : :
470 : : // Categorize mesh data (tets, node coordinates, and boundary data) by target
471 : : // chare based on which chare the partitioner assigned elements (tets) to
472 : : std::unordered_map< int, MeshData > chmesh;
473 [ + + ]: 310333 : for (std::size_t e=0; e<target.size(); ++e) {
474 : : // Construct a tetrahedron with global node ids
475 [ + - ]: 309806 : tk::UnsMesh::Tet t{{ m_ginpoel[e*4+0], m_ginpoel[e*4+1],
476 : 309806 : m_ginpoel[e*4+2], m_ginpoel[e*4+3] }};
477 : : // Categorize tetrahedron (domain element) connectivity
478 [ + - ][ + - ]: 309806 : auto& mesh = chmesh[ static_cast<int>(target[e]) ];
479 : : auto& inpoel = std::get< 0 >( mesh );
480 : : inpoel.insert( end(inpoel), begin(t), end(t) );
481 : : // Categorize mesh block ids and associated local tet ids (local tet ids
482 : : // basically correspond to the inpoel entries updated in line above).
483 : : [[maybe_unused]] bool added(false);
484 : : auto& elblock = std::get< 3 >( mesh );
485 [ + + ]: 627404 : for (const auto& [blid, elset] : m_elemBlockId) {
486 [ + + ]: 317598 : if (elset.find(e) != elset.end()) {
487 [ + - ]: 309806 : elblock.push_back(blid);
488 : : added = true;
489 : : }
490 : : }
491 : : Assert(added, "Tet element not found in any block after partitioning");
492 : : // Categorize boundary face connectivity
493 : : auto& bconn = std::get< 1 >( mesh );
494 : : std::array<Face,4> face{{ {{t[0],t[2],t[1]}}, {{t[0],t[1],t[3]}},
495 : 309806 : {{t[0],t[3],t[2]}}, {{t[1],t[2],t[3]}} }};
496 [ + + ]: 1549030 : for (const auto& f : face) {
497 : : auto it = faceside.find( f );
498 [ + + ]: 1239224 : if (it != end(faceside)) {
499 [ + - ]: 134482 : auto& s = bconn[ it->second ];
500 : : s.insert( end(s), begin(f), end(f) );
501 : : }
502 : : }
503 : : // Categorize boundary node lists
504 : : auto& bnode = std::get< 2 >( mesh );
505 [ + + ]: 1549030 : for (const auto& n : t) {
506 : : auto it = nodeside.find( n );
507 [ + + ]: 1239224 : if (it != end(nodeside))
508 [ + + ][ + - ]: 641318 : for (auto s : it->second)
509 [ + - ]: 335032 : bnode[ s ].push_back( n );
510 : : }
511 : : }
512 : :
513 : : // Make boundary node lists unique per side set
514 [ + + ]: 7672 : for (auto& c : chmesh)
515 [ + + ]: 13648 : for (auto& n : std::get<2>(c.second))
516 [ + - ]: 6503 : tk::unique( n.second );
517 : :
518 : : // Make sure all compute nodes have target chares assigned
519 : : Assert( !chmesh.empty(), "No elements have been assigned to a chare" );
520 : :
521 : : // This check should always be done, hence ErrChk and not Assert, as it
522 : : // can result from particular pathological combinations of (1) too large
523 : : // degree of virtualization, (2) too many compute nodes, and/or (3) too small
524 : : // of a mesh and not due to programmer error.
525 [ + + ]: 7672 : for(const auto& c : chmesh)
526 [ - + ][ - - ]: 7145 : ErrChk( !std::get<0>(c.second).empty(),
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
527 : : "Overdecomposition of the mesh is too large compared to the "
528 : : "number of work units computed based on the degree of "
529 : : "virtualization desired. As a result, there would be at least "
530 : : "one work unit with no mesh elements to work on, i.e., nothing "
531 : : "to do. Solution 1: decrease the virtualization to a lower "
532 : : "value using the command-line argument '-u'. Solution 2: "
533 : : "decrease the number processing elements (PEs and/or compute "
534 : : "nodes) using the charmrun command-line argument '+pN' where N is "
535 : : "the number of PEs (or in SMP-mode in combination with +ppn to "
536 : : "reduce the number of compute nodes), which implicitly increases "
537 : : "the size (and thus decreases the number) of work units.)" );
538 : :
539 : 527 : return chmesh;
540 : : }
541 : :
542 : : tk::UnsMesh::CoordMap
543 [ + - ]: 7145 : Partitioner::coordmap( const std::vector< std::size_t >& inpoel )
544 : : // *****************************************************************************
545 : : // Extract coordinates associated to global nodes of a mesh chunk
546 : : //! \param[in] inpoel Mesh connectivity
547 : : //! \return Map storing the coordinates of unique nodes associated to global
548 : : //! node IDs in mesh given by inpoel
549 : : // *****************************************************************************
550 : : {
551 : : Assert( inpoel.size() % 4 == 0, "Incomplete mesh connectivity" );
552 : :
553 : : tk::UnsMesh::CoordMap map;
554 : :
555 [ + - ][ + + ]: 266472 : for (auto g : tk::uniquecopy(inpoel)) {
556 [ + - ]: 252182 : auto i = tk::cref_find( m_lid, g );
557 : : auto& c = map[g];
558 : 252182 : c[0] = m_coord[0][i];
559 : 252182 : c[1] = m_coord[1][i];
560 : 252182 : c[2] = m_coord[2][i];
561 : : }
562 : :
563 : : Assert( tk::uniquecopy(inpoel).size() == map.size(), "Size mismatch" );
564 : :
565 : 7145 : return map;
566 : : }
567 : :
568 : : void
569 : 527 : Partitioner::distribute( std::unordered_map< int, MeshData >&& mesh )
570 : : // *****************************************************************************
571 : : // Distribute mesh to target compute nodes after mesh partitioning
572 : : //! \param[in] mesh Mesh data categorized by target by target chares
573 : : // *****************************************************************************
574 : : {
575 : 527 : auto dist = distribution( m_nchare );
576 : :
577 : : // Extract mesh data whose chares are on ("owned by") this compute node
578 [ + + ]: 2235 : for (int c=0; c<dist[1]; ++c) {
579 : 1708 : auto chid = CkMyNode() * dist[0] + c; // compute owned chare ID
580 : : const auto it = mesh.find( chid ); // attempt to find its mesh data
581 [ + + ]: 1708 : if (it != end(mesh)) { // if found
582 : : // Store own tetrahedron connectivity
583 : : const auto& inpoel = std::get<0>( it->second );
584 : : auto& inp = m_chinpoel[ chid ]; // will store own mesh connectivity
585 : 1495 : inp.insert( end(inp), begin(inpoel), end(inpoel) );
586 : : // Store own mesh block id and associated tet ids
587 : : const auto& elblock = std::get<3>( it->second );
588 : : auto& cheb = m_chelemblockid[ chid ];
589 : 1495 : cheb.insert( end(cheb), begin(elblock), end(elblock) );
590 : : // Store own node coordinates
591 : : auto& chcm = m_chcoordmap[ chid ]; // will store own node coordinates
592 : 1495 : auto cm = coordmap( inpoel ); // extract node coordinates
593 : : chcm.insert( begin(cm), end(cm) ); // concatenate node coords
594 : : // Store own boundary face connectivity
595 : : const auto& bconn = std::get<1>( it->second );
596 : : auto& bface = m_chbface[ chid ]; // will store own boundary faces
597 [ + - ]: 1495 : auto& t = m_chtriinpoel[ chid ]; // wil store own boundary face conn
598 : : auto& f = m_nface[ chid ]; // use counter for chare
599 [ + + ]: 4339 : for (const auto& [ setid, faceids ] : bconn) {
600 [ + - ]: 2844 : auto& b = bface[ setid ];
601 [ + + ]: 80384 : for (std::size_t i=0; i<faceids.size()/3; ++i) {
602 [ + - ]: 77540 : b.push_back( f++ );
603 [ + - ]: 77540 : t.push_back( faceids[i*3+0] );
604 [ + - ]: 77540 : t.push_back( faceids[i*3+1] );
605 [ + - ]: 77540 : t.push_back( faceids[i*3+2] );
606 : : }
607 : : }
608 : : // Store own boundary node lists
609 : : const auto& bnode = std::get<2>( it->second );
610 : : auto& nodes = m_chbnode[ chid ]; // will store own boundary nodes
611 [ + + ]: 2988 : for (const auto& [ setid, nodeids ] : bnode) {
612 [ + - ]: 1493 : auto& b = nodes[ setid ];
613 [ + - ]: 1493 : b.insert( end(b), begin(nodeids), end(nodeids) );
614 : : }
615 : : // Remove chare ID and mesh data
616 : : mesh.erase( it );
617 : : }
618 : : Assert( mesh.find(chid) == end(mesh), "Not all owned mesh data stored" );
619 : : }
620 : :
621 : : // Construct export map (associating mesh connectivities with global node
622 : : // indices and node coordinates) for mesh chunks associated to chare IDs
623 : : // owned by chares we do not own.
624 : : std::unordered_map< int, // target compute node
625 : : std::unordered_map< int, // chare ID
626 : : std::tuple<
627 : : // (domain-element) tetrahedron connectivity
628 : : std::vector< std::size_t >,
629 : : // (domain) node IDs & coordinates
630 : : tk::UnsMesh::CoordMap,
631 : : // boundary side set + face connectivity
632 : : std::unordered_map< int, std::vector< std::size_t > >,
633 : : // boundary side set + node list
634 : : std::unordered_map< int, std::vector< std::size_t > >,
635 : : // Mesh block ids + local tet ids
636 : : std::vector< std::size_t >
637 : : > > > exp;
638 : :
639 [ + + ]: 6177 : for (const auto& c : mesh)
640 [ + - ][ + - ]: 11300 : exp[ node(c.first) ][ c.first ] =
[ + - ]
641 : 5650 : std::make_tuple( std::get<0>(c.second),
642 [ + - ]: 11300 : coordmap(std::get<0>(c.second)),
643 : : std::get<1>(c.second),
644 : : std::get<2>(c.second),
645 : : std::get<3>(c.second) );
646 : :
647 : : // Export chare IDs and mesh we do not own to fellow compute nodes
648 [ + + ]: 527 : if (exp.empty()) {
649 [ - + ][ - - ]: 77 : if ( g_inputdeck.get< tag::cmd, tag::feedback >() ) m_host.pedistributed();
650 [ + - ]: 77 : contribute( sizeof(std::size_t), &m_meshid, CkReduction::nop,
651 : : m_cbp.get< tag::distributed >() );
652 : : } else {
653 : 450 : m_ndist += exp.size();
654 [ + + ]: 2442 : for (const auto& [ targetnode, chunk ] : exp)
655 [ + - ][ + - ]: 5976 : thisProxy[ targetnode ].addMesh( CkMyNode(), chunk );
656 : : }
657 : 527 : }
658 : :
659 : : std::array< int, 2 >
660 : 1054 : Partitioner::distribution( int npart ) const
661 : : // *****************************************************************************
662 : : // Compute chare (partition) distribution
663 : : //! \param[in] npart Total number of chares (partitions) to distribute
664 : : //! \return Chunksize, i.e., number of chares per all compute nodes except the
665 : : //! last one, and the number of chares for this compute node.
666 : : //! \details Chare ids are distributed to compute nodes in a linear continguous
667 : : //! order with the last compute node taking the remainder if the number of
668 : : //! compute nodes is not divisible by the number chares. For example, if
669 : : //! nchare=7 and nnode=3, the chare distribution is node0: 0 1, node1: 2 3,
670 : : //! and node2: 4 5 6. As a result of this distribution, all compute nodes will
671 : : //! have their chare-categorized element connectivity filled with the global
672 : : //! mesh node IDs associated to the Charm++ chare IDs each compute node owns.
673 : : // *****************************************************************************
674 : : {
675 : 1054 : auto chunksize = npart / CkNumNodes();
676 : : auto mynchare = chunksize;
677 [ + + ]: 1054 : if (CkMyNode() == CkNumNodes()-1) mynchare += npart % CkNumNodes();
678 : 1054 : return {{ chunksize, mynchare }};
679 : : }
680 : :
681 [ - - ]: 0 : Partitioner::~Partitioner()
682 : : // *****************************************************************************
683 : : // Destructor
684 : : // *****************************************************************************
685 : : {
686 : : // The following call has to be made on all MPI ranks to free all resources.
687 : : // see https://github.com/trilinos/Trilinos/issues/11197#issuecomment-1301325163
688 [ - - ]: 0 : if (m_meshid == 0) Kokkos::finalize();
689 : 0 : }
690 : :
691 : : #include "NoWarning/partitioner.def.h"
|
