Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/Refiner.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 Mesh refiner for interfacing the mesh refinement library
9 : : \see Refiner.h for more info.
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include <vector>
14 : : #include <algorithm>
15 : :
16 : : #include "Refiner.hpp"
17 : : #include "Reorder.hpp"
18 : : #include "AMR/mesh_adapter.hpp"
19 : : #include "AMR/Error.hpp"
20 : : #include "Inciter/InputDeck/InputDeck.hpp"
21 : : #include "CGPDE.hpp"
22 : : #include "DGPDE.hpp"
23 : : #include "FVPDE.hpp"
24 : : #include "DerivedData.hpp"
25 : : #include "UnsMesh.hpp"
26 : : #include "Centering.hpp"
27 : : #include "Around.hpp"
28 : : #include "Sorter.hpp"
29 : : #include "Discretization.hpp"
30 : :
31 : : namespace inciter {
32 : :
33 : : extern ctr::InputDeck g_inputdeck;
34 : : extern ctr::InputDeck g_inputdeck_defaults;
35 : : extern std::vector< CGPDE > g_cgpde;
36 : : extern std::vector< DGPDE > g_dgpde;
37 : : extern std::vector< FVPDE > g_fvpde;
38 : :
39 : : } // inciter::
40 : :
41 : : using inciter::Refiner;
42 : :
43 : 1943 : Refiner::Refiner( std::size_t meshid,
44 : : const CProxy_Transporter& transporter,
45 : : const CProxy_Sorter& sorter,
46 : : const tk::CProxy_MeshWriter& meshwriter,
47 : : const std::vector< Scheme >& scheme,
48 : : const tk::RefinerCallback& cbr,
49 : : const tk::SorterCallback& cbs,
50 : : const std::vector< std::size_t >& ginpoel,
51 : : const tk::UnsMesh::CoordMap& coordmap,
52 : : const std::map< int, std::vector< std::size_t > >& bface,
53 : : const std::vector< std::size_t >& triinpoel,
54 : : const std::map< int, std::vector< std::size_t > >& bnode,
55 : : const std::vector< std::size_t >& elemblid,
56 : 1943 : int nchare ) :
57 : : m_meshid( meshid ),
58 : : m_ncit(0),
59 : : m_host( transporter ),
60 : : m_sorter( sorter ),
61 : : m_meshwriter( meshwriter ),
62 : : m_scheme( scheme ),
63 : : m_cbr( cbr ),
64 : : m_cbs( cbs ),
65 : : m_ginpoel( ginpoel ),
66 : : m_el( tk::global2local( ginpoel ) ), // fills m_inpoel, m_gid, m_lid
67 : : m_coordmap( coordmap ),
68 : : m_coord( flatcoord(coordmap) ),
69 : : m_bface( bface ),
70 : : m_bnode( bnode ),
71 : : m_triinpoel( triinpoel ),
72 : : m_elemblockid(),
73 : : m_nchare( nchare ),
74 : : m_mode( RefMode::T0REF ),
75 : : m_initref( g_inputdeck.get< tag::amr, tag::initial >() ),
76 [ + - ]: 1943 : m_ninitref( g_inputdeck.get< tag::amr, tag::initial >().size() ),
77 : 1943 : m_refiner( g_inputdeck.get< tag::amr, tag::maxlevels >(), m_inpoel ),
78 : : m_nref( 0 ),
79 : : m_nbnd( 0 ),
80 : : m_extra( 0 ),
81 : : m_ch(),
82 : : m_edgech(),
83 : : m_chedge(),
84 : : m_localEdgeData(),
85 : : m_remoteEdgeData(),
86 : : m_nodeCommMap(),
87 : : m_oldTets(),
88 : : m_addedNodes(),
89 : : m_addedTets(),
90 : : m_removedNodes(),
91 : : m_amrNodeMap(),
92 : : m_oldntets( 0 ),
93 : : m_coarseBndFaces(),
94 : : m_coarseBndNodes(),
95 : : m_coarseBlkElems(),
96 : : m_rid( m_coord[0].size() ),
97 : : // m_oldrid(),
98 : : m_lref( m_rid.size() ),
99 : : // m_oldparent(),
100 : : m_writeCallback(),
101 : : m_outref_ginpoel(),
102 : : m_outref_el(),
103 : : m_outref_coord(),
104 : : m_outref_addedNodes(),
105 : : m_outref_addedTets(),
106 : : m_outref_nodeCommMap(),
107 : : m_outref_bface(),
108 : : m_outref_bnode(),
109 [ + - ][ + - ]: 7772 : m_outref_triinpoel()
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ]
110 : : // *****************************************************************************
111 : : // Constructor
112 : : //! \param[in] meshid Mesh ID
113 : : //! \param[in] transporter Transporter (host) proxy
114 : : //! \param[in] sorter Mesh reordering (sorter) proxy
115 : : //! \param[in] meshwriter Mesh writer proxy
116 : : //! \param[in] scheme Discretization schemes (one per mesh)
117 : : //! \param[in] cbr Charm++ callbacks for Refiner
118 : : //! \param[in] cbs Charm++ callbacks for Sorter
119 : : //! \param[in] ginpoel Mesh connectivity (this chare) using global node IDs
120 : : //! \param[in] coordmap Mesh node coordinates (this chare) for global node IDs
121 : : //! \param[in] bface File-internal elem ids of side sets
122 : : //! \param[in] triinpoel Triangle face connectivity with global IDs
123 : : //! \param[in] bnode Node lists of side sets
124 : : //! \param[in] elemblid Mesh block ids associated to local tet ids
125 : : //! \param[in] nchare Total number of refiner chares (chare array elements)
126 : : // *****************************************************************************
127 : : {
128 : : Assert( !m_ginpoel.empty(), "No elements assigned to refiner chare" );
129 : : Assert( tk::positiveJacobians( m_inpoel, m_coord ),
130 : : "Input mesh to Refiner Jacobian non-positive" );
131 : : Assert( !tk::leakyPartition(
132 : : tk::genEsuelTet( m_inpoel, tk::genEsup(m_inpoel,4) ),
133 : : m_inpoel, m_coord ),
134 : : "Input mesh to Refiner leaky" );
135 : :
136 : : // Construct data structure assigning sets of element ids to mesh blocks
137 [ + + ]: 448500 : for (std::size_t ie=0; ie<elemblid.size(); ++ie) {
138 : : m_elemblockid[elemblid[ie]].insert(ie);
139 : : }
140 : :
141 : : #if not defined(__INTEL_COMPILER) || defined(NDEBUG)
142 : : // The above ifdef skips running the conformity test with the intel compiler
143 : : // in debug mode only. This is necessary because in tk::conforming(), filling
144 : : // up the map can fail with some meshes (only in parallel), e.g., tube.exo,
145 : : // used by some regression tests, due to the intel compiler generating some
146 : : // garbage incorrect code - only in debug, only in parallel, only with that
147 : : // mesh.
148 : : Assert( tk::conforming( m_inpoel, m_coord, true, m_rid ),
149 : : "Input mesh to Refiner not conforming" );
150 : : #endif
151 : :
152 : : // Generate local -> refiner lib node id map and its inverse
153 [ + - ]: 1943 : libmap();
154 : :
155 : : // Reverse initial mesh refinement type list (will pop from back)
156 : : std::reverse( begin(m_initref), end(m_initref) );
157 : :
158 : : // Generate boundary data structures for coarse mesh
159 [ + - ]: 1943 : coarseMesh();
160 : :
161 : : // If initial mesh refinement is configured, start initial mesh refinement.
162 : : // See also tk::grm::check_amr_errors in Control/Inciter/InputDeck/Ggrammar.h.
163 [ + + ][ + - ]: 1943 : if (g_inputdeck.get< tag::amr, tag::t0ref >() && m_ninitref > 0)
164 [ + - ]: 38 : t0ref();
165 : : else
166 [ + - ]: 1905 : endt0ref();
167 : 1943 : }
168 : :
169 : : void
170 : 1943 : Refiner::libmap()
171 : : // *****************************************************************************
172 : : // (Re-)generate local -> refiner lib node id map and its inverse
173 : : // *****************************************************************************
174 : : {
175 : : // Fill initial (matching) mapping between local and refiner node ids
176 : : std::iota( begin(m_rid), end(m_rid), 0 );
177 : :
178 : : // Fill in inverse, mapping refiner to local node ids
179 : : std::size_t i = 0;
180 [ + + ]: 164234 : for (auto r : m_rid) m_lref[r] = i++;
181 : 1943 : }
182 : :
183 : : void
184 : 1961 : Refiner::coarseMesh()
185 : : // *****************************************************************************
186 : : // (Re-)generate side set and block data structures for coarse mesh
187 : : // *****************************************************************************
188 : : {
189 : : // Generate unique set of faces for each side set of the input (coarsest) mesh
190 : : m_coarseBndFaces.clear();
191 [ + + ]: 4820 : for (const auto& [ setid, faceids ] : m_bface) {
192 : : auto& faces = m_coarseBndFaces[ setid ];
193 [ + + ]: 169839 : for (auto f : faceids) {
194 : 166980 : faces.insert(
195 : 166980 : {{{ m_triinpoel[f*3+0], m_triinpoel[f*3+1], m_triinpoel[f*3+2] }}} );
196 : : }
197 : : }
198 : :
199 : : // Generate unique set of nodes for each side set of the input (coarsest) mesh
200 : : m_coarseBndNodes.clear();
201 [ + + ]: 2499 : for (const auto& [ setid, nodes ] : m_bnode) {
202 : : m_coarseBndNodes[ setid ].insert( begin(nodes), end(nodes) );
203 : : }
204 : :
205 : : // Generate set of elements for each mesh block of the input (coarsest) mesh
206 : : m_coarseBlkElems.clear();
207 [ + + ]: 3940 : for (const auto& [blid, elids] : m_elemblockid) {
208 [ + + ][ + - ]: 453034 : for (auto ie : elids) {
209 [ + - ]: 451055 : m_coarseBlkElems[blid].insert( {{{m_inpoel[ie*4+0], m_inpoel[ie*4+1],
210 [ + - ]: 451055 : m_inpoel[ie*4+2], m_inpoel[ie*4+3]}}} );
211 : : }
212 : : }
213 : 1961 : }
214 : :
215 : : void
216 : 1943 : Refiner::sendProxy()
217 : : // *****************************************************************************
218 : : // Send Refiner proxy to Discretization objects
219 : : //! \details This should be called when bound Discretization chare array
220 : : //! elements have already been created.
221 : : // *****************************************************************************
222 : : {
223 : : // Make sure (bound) Discretization chare is already created and accessible
224 : : Assert( m_scheme[m_meshid].disc()[thisIndex].ckLocal() != nullptr,
225 : : "About to dereference nullptr" );
226 : :
227 : : // Pass Refiner Charm++ chare proxy to fellow (bound) Discretization object
228 [ + - ]: 3886 : m_scheme[m_meshid].disc()[thisIndex].ckLocal()->setRefiner( thisProxy );
229 : 1943 : }
230 : :
231 : : void
232 : 18 : Refiner::reorder()
233 : : // *****************************************************************************
234 : : // Query Sorter and update local mesh with the reordered one
235 : : // *****************************************************************************
236 : : {
237 : 18 : m_sorter[thisIndex].ckLocal()->
238 [ + - ]: 36 : mesh( m_ginpoel, m_coordmap, m_triinpoel, m_bnode );
239 : :
240 : : // Update local mesh data based on data just received from Sorter
241 : 18 : m_el = tk::global2local( m_ginpoel ); // fills m_inpoel, m_gid, m_lid
242 : 18 : m_coord = flatcoord( m_coordmap );
243 : :
244 : : // Re-generate boundary data structures for coarse mesh
245 : 18 : coarseMesh();
246 : :
247 : : // WARNING: This re-creates the AMR lib which is probably not what we
248 : : // ultimately want, beacuse this deletes its history recorded during initial
249 : : // (t<0) refinement. However, this appears to correctly update the local mesh
250 : : // based on the reordered one (from Sorter) at least when t0ref is off.
251 : 18 : m_refiner = AMR::mesh_adapter_t(
252 : 18 : g_inputdeck.get< tag::amr, tag::maxlevels >(), m_inpoel );
253 : 18 : }
254 : :
255 : : tk::UnsMesh::Coords
256 : 1961 : Refiner::flatcoord( const tk::UnsMesh::CoordMap& coordmap )
257 : : // *****************************************************************************
258 : : // Generate flat coordinate data from coordinate map
259 : : //! \param[in] coordmap Coordinates associated to global node IDs of mesh chunk
260 : : //! \return Flat coordinate data
261 : : // *****************************************************************************
262 : : {
263 : : tk::UnsMesh::Coords coord;
264 : :
265 : : // Convert node coordinates associated to global node IDs to a flat vector
266 : : auto npoin = coordmap.size();
267 : : Assert( m_gid.size() == npoin, "Size mismatch" );
268 [ + - ]: 1961 : coord[0].resize( npoin );
269 [ + - ]: 1961 : coord[1].resize( npoin );
270 [ + - ]: 1961 : coord[2].resize( npoin );
271 [ + + ]: 176641 : for (const auto& [ gid, coords ] : coordmap) {
272 : 174680 : auto i = tk::cref_find( m_lid, gid );
273 : : Assert( i < npoin, "Indexing out of coordinate map" );
274 : 174680 : coord[0][i] = coords[0];
275 : 174680 : coord[1][i] = coords[1];
276 : 174680 : coord[2][i] = coords[2];
277 : : }
278 : :
279 : 1961 : return coord;
280 : : }
281 : :
282 : : void
283 : 0 : Refiner::dtref( const std::map< int, std::vector< std::size_t > >& bface,
284 : : const std::map< int, std::vector< std::size_t > >& bnode,
285 : : const std::vector< std::size_t >& triinpoel )
286 : : // *****************************************************************************
287 : : // Start mesh refinement (during time stepping, t>0)
288 : : //! \param[in] bface Boundary-faces mapped to side set ids
289 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
290 : : //! \param[in] triinpoel Boundary-face connectivity
291 : : // *****************************************************************************
292 : : {
293 : 0 : m_mode = RefMode::DTREF;
294 : :
295 : : // Update boundary node lists
296 : : m_bface = bface;
297 : : m_bnode = bnode;
298 : 0 : m_triinpoel = tk::remap(triinpoel, m_gid);
299 : :
300 : 0 : start();
301 : 0 : }
302 : :
303 : : void
304 : 66 : Refiner::outref( const std::map< int, std::vector< std::size_t > >& bface,
305 : : const std::map< int, std::vector< std::size_t > >& bnode,
306 : : const std::vector< std::size_t >& triinpoel,
307 : : CkCallback c,
308 : : RefMode mode )
309 : : // *****************************************************************************
310 : : // Start mesh refinement (for field output)
311 : : //! \param[in] bface Boundary-faces mapped to side set ids
312 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
313 : : //! \param[in] triinpoel Boundary-face connectivity
314 : : //! \param[in] c Function to continue with after the writing field output
315 : : //! \param[in] mode Refinement mode
316 : : // *****************************************************************************
317 : : {
318 : 66 : m_mode = mode;
319 : :
320 : 66 : m_writeCallback = c;
321 : :
322 : : // Update boundary node lists
323 : : m_bface = bface;
324 : : m_bnode = bnode;
325 : 66 : m_triinpoel = triinpoel;
326 : :
327 : 66 : start();
328 : 66 : }
329 : :
330 : : void
331 : 111 : Refiner::t0ref()
332 : : // *****************************************************************************
333 : : // Output mesh to file before a new step mesh refinement
334 : : // *****************************************************************************
335 : : {
336 : : Assert( m_ninitref > 0, "No initial mesh refinement steps configured" );
337 : : // Output initial mesh to file
338 [ + + ]: 111 : auto l = m_ninitref - m_initref.size(); // num initref steps completed
339 : 111 : auto t0 = g_inputdeck.get< tag::t0 >();
340 [ + + ]: 111 : if (l == 0) {
341 [ + - ][ + - ]: 76 : writeMesh( "t0ref", l, t0-1.0,
342 [ + - ][ - + ]: 114 : CkCallback( CkIndex_Refiner::start(), thisProxy[thisIndex] ) );
[ - - ]
343 : : } else {
344 : 73 : start();
345 : : }
346 : 111 : }
347 : :
348 : : void
349 : 177 : Refiner::start()
350 : : // *****************************************************************************
351 : : // Start new step of mesh refinement
352 : : // *****************************************************************************
353 : : {
354 : 177 : m_extra = 0;
355 : : m_ch.clear();
356 : : m_remoteEdgeData.clear();
357 : : m_remoteEdges.clear();
358 : :
359 : 177 : updateEdgeData();
360 : :
361 : : // Generate and communicate boundary edges
362 : 177 : bndEdges();
363 : 177 : }
364 : :
365 : : void
366 : 177 : Refiner::bndEdges()
367 : : // *****************************************************************************
368 : : // Generate boundary edges and send them to all chares
369 : : //! \details Extract edges on the boundary only. The boundary edges (shared by
370 : : //! multiple chares) will be agreed on a refinement that yields a conforming
371 : : //! mesh across chares boundaries.
372 : : // *****************************************************************************
373 : : {
374 : : // Compute the number of edges (chunksize) a chare will respond to when
375 : : // computing shared edges
376 : 177 : auto N = static_cast< std::size_t >( m_nchare );
377 [ + - ]: 177 : std::size_t chunksize = std::numeric_limits< std::size_t >::max() / N;
378 : :
379 : : // Generate boundary edges of our mesh chunk
380 : : std::unordered_map< int, EdgeSet > chbedges;
381 [ + - ]: 354 : auto esup = tk::genEsup( m_inpoel, 4 ); // elements surrounding points
382 [ + - ]: 177 : auto esuel = tk::genEsuelTet( m_inpoel, esup ); // elems surrounding elements
383 [ + + ]: 294172 : for (std::size_t e=0; e<esuel.size()/4; ++e) {
384 : 293995 : auto mark = e*4;
385 [ + + ]: 1469975 : for (std::size_t f=0; f<4; ++f) {
386 [ + + ]: 1175980 : if (esuel[mark+f] == -1) {
387 [ + - ]: 111628 : auto A = m_ginpoel[ mark+tk::lpofa[f][0] ];
388 : 111628 : auto B = m_ginpoel[ mark+tk::lpofa[f][1] ];
389 : 111628 : auto C = m_ginpoel[ mark+tk::lpofa[f][2] ];
390 : : Assert( m_lid.find( A ) != end(m_lid), "Local node ID not found" );
391 : : Assert( m_lid.find( B ) != end(m_lid), "Local node ID not found" );
392 : : Assert( m_lid.find( C ) != end(m_lid), "Local node ID not found" );
393 : : // assign edges to bins a single chare will respond to when computing
394 : : // shared edges
395 : 111628 : auto bin = A / chunksize;
396 : : Assert( bin < N, "Will index out of number of chares" );
397 [ + - ][ + - ]: 111628 : chbedges[ static_cast<int>(bin) ].insert( {A,B} );
398 : 111628 : bin = B / chunksize;
399 : : Assert( bin < N, "Will index out of number of chares" );
400 [ + - ][ + - ]: 111628 : chbedges[ static_cast<int>(bin) ].insert( {B,C} );
401 : 111628 : bin = C / chunksize;
402 : : Assert( bin < N, "Will index out of number of chares" );
403 [ + - ][ + - ]: 111628 : chbedges[ static_cast<int>(bin) ].insert( {C,A} );
404 : : }
405 : : }
406 : : }
407 : :
408 : : // Send edges in bins to chares that will compute shared edges
409 : 177 : m_nbnd = chbedges.size();
410 [ - + ]: 177 : if (m_nbnd == 0)
411 [ - - ]: 0 : contribute( sizeof(std::size_t), &m_meshid, CkReduction::nop,
412 : : m_cbr.get< tag::queried >() );
413 : : else
414 [ + + ]: 488 : for (const auto& [ targetchare, bndedges ] : chbedges)
415 [ + - ][ + - ]: 622 : thisProxy[ targetchare ].query( thisIndex, bndedges );
[ - - ]
416 : 177 : }
417 : :
418 : : void
419 : 311 : Refiner::query( int fromch, const EdgeSet& edges )
420 : : // *****************************************************************************
421 : : // Incoming query for a list boundary edges for which this chare compiles
422 : : // shared edges
423 : : //! \param[in] fromch Sender chare ID
424 : : //! \param[in] edges Chare-boundary edge list from another chare
425 : : // *****************************************************************************
426 : : {
427 : : // Store incoming edges in edge->chare and its inverse, chare->edge, maps
428 [ + + ]: 403949 : for (const auto& e : edges) m_edgech[ e ].push_back( fromch );
429 : : m_chedge[ fromch ].insert( begin(edges), end(edges) );
430 : : // Report back to chare message received from
431 [ + - ]: 311 : thisProxy[ fromch ].recvquery();
432 : 311 : }
433 : :
434 : : void
435 : 311 : Refiner::recvquery()
436 : : // *****************************************************************************
437 : : // Receive receipt of boundary edge lists to query
438 : : // *****************************************************************************
439 : : {
440 [ + + ]: 311 : if (--m_nbnd == 0)
441 : 177 : contribute( sizeof(std::size_t), &m_meshid, CkReduction::nop,
442 : : m_cbr.get< tag::queried >() );
443 : 311 : }
444 : :
445 : : void
446 : 177 : Refiner::response()
447 : : // *****************************************************************************
448 : : // Respond to boundary edge list queries
449 : : // *****************************************************************************
450 : : {
451 : : std::unordered_map< int, std::vector< int > > exp;
452 : :
453 : : // Compute shared edges whose chare ids will be sent back to querying chares
454 [ + + ]: 488 : for (const auto& [ neighborchare, bndedges ] : m_chedge) {
455 : : auto& e = exp[ neighborchare ];
456 [ + + ]: 202130 : for (const auto& ed : bndedges)
457 [ + + ]: 429236 : for (auto d : tk::cref_find(m_edgech,ed))
458 [ + + ]: 227417 : if (d != neighborchare)
459 [ + - ]: 25598 : e.push_back( d );
460 : : }
461 : :
462 : : // Send chare ids of shared edges to chares that issued a query to us. Shared
463 : : // boundary edges assigned to chare ids sharing the boundary edge were
464 : : // computed above for those chares that queried this map from us. These
465 : : // boundary edges form a distributed table and we only work on a chunk of it.
466 : : // Note that we only send data back to those chares that have queried us. The
467 : : // receiving sides do not know in advance if they receive messages or not.
468 : : // Completion is detected by having the receiver respond back and counting
469 : : // the responses on the sender side, i.e., this chare.
470 : 177 : m_nbnd = exp.size();
471 [ + + ]: 177 : if (m_nbnd == 0)
472 [ + - ]: 50 : contribute( sizeof(std::size_t), &m_meshid, CkReduction::nop,
473 : : m_cbr.get< tag::responded >() );
474 : : else
475 [ + + ]: 438 : for (const auto& [ targetchare, bndedges ] : exp)
476 [ + - ][ + - ]: 622 : thisProxy[ targetchare ].bnd( thisIndex, bndedges );
477 : 177 : }
478 : :
479 : : void
480 : 311 : Refiner::bnd( int fromch, const std::vector< int >& chares )
481 : : // *****************************************************************************
482 : : // Receive shared boundary edges for our mesh chunk
483 : : //! \param[in] fromch Sender chare ID
484 : : //! \param[in] chares Chare ids we share edges with
485 : : // *****************************************************************************
486 : : {
487 : : // Store chare ids we share edges with
488 : : m_ch.insert( begin(chares), end(chares) );
489 : :
490 : : // Report back to chare message received from
491 [ + - ]: 311 : thisProxy[ fromch ].recvbnd();
492 : 311 : }
493 : :
494 : : void
495 : 311 : Refiner::recvbnd()
496 : : // *****************************************************************************
497 : : // Receive receipt of shared boundary edges
498 : : // *****************************************************************************
499 : : {
500 [ + + ]: 311 : if (--m_nbnd == 0)
501 : 127 : contribute( sizeof(std::size_t), &m_meshid, CkReduction::nop,
502 : : m_cbr.get< tag::responded >() );
503 : 311 : }
504 : :
505 : : void
506 : 177 : Refiner::refine()
507 : : // *****************************************************************************
508 : : // Do a single step of mesh refinement (really, only tag edges)
509 : : //! \details During initial (t<0, t0ref) mesh refinement, this is a single step
510 : : //! in a potentially multiple-entry list of initial adaptive mesh refinement
511 : : //! steps. Distribution of the chare-boundary edges must have preceded this
512 : : //! step, so that boundary edges (shared by multiple chares) can agree on a
513 : : //! refinement that yields a conforming mesh across chare boundaries.
514 : : //!
515 : : //! During-timestepping (t>0, dtref) mesh refinement this is called once, as
516 : : //! we only do a single step during time stepping.
517 : : //!
518 : : //! During field-output (outref) mesh refinement, this may be called multiple
519 : : //! times, depending the number of refinement levels needed for the field
520 : : //! output.
521 : : // *****************************************************************************
522 : : {
523 : : // Free memory used for computing shared boundary edges
524 : 177 : tk::destroy( m_edgech );
525 : 177 : tk::destroy( m_chedge );
526 : :
527 : : // Perform leak test on old mesh
528 : : Assert( !tk::leakyPartition(
529 : : tk::genEsuelTet( m_inpoel, tk::genEsup(m_inpoel,4) ),
530 : : m_inpoel, m_coord ),
531 : : "Mesh partition before refinement leaky" );
532 : :
533 [ + + ]: 177 : if (m_mode == RefMode::T0REF) {
534 : :
535 : : // Refine mesh based on next initial refinement type
536 [ + - ]: 111 : if (!m_initref.empty()) {
537 : 111 : auto r = m_initref.back(); // consume (reversed) list from its back
538 [ + + ]: 111 : if (r == ctr::AMRInitialType::UNIFORM)
539 : 57 : uniformRefine();
540 [ + + ]: 54 : else if (r == ctr::AMRInitialType::UNIFORM_DEREFINE)
541 : 36 : uniformDeRefine();
542 [ + + ]: 18 : else if (r == ctr::AMRInitialType::INITIAL_CONDITIONS)
543 : 16 : errorRefine();
544 [ + - ]: 2 : else if (r == ctr::AMRInitialType::COORDINATES)
545 : 2 : coordRefine();
546 [ - - ]: 0 : else if (r == ctr::AMRInitialType::EDGELIST)
547 : 0 : edgelistRefine();
548 [ - - ][ - - ]: 0 : else Throw( "Initial AMR type not implemented" );
[ - - ][ - - ]
[ - - ][ - - ]
549 : : }
550 : :
551 [ - + ]: 66 : } else if (m_mode == RefMode::DTREF) {
552 : :
553 [ - - ]: 0 : if (g_inputdeck.get< tag::amr, tag::dtref_uniform >())
554 : 0 : uniformRefine();
555 : : else
556 : 0 : errorRefine();
557 : :
558 [ + + ]: 66 : } else if (m_mode == RefMode::OUTREF) {
559 : :
560 : 33 : uniformRefine();
561 : :
562 [ + - ]: 33 : } else if (m_mode == RefMode::OUTDEREF) {
563 : :
564 : 33 : uniformDeRefine();
565 : :
566 [ - - ][ - - ]: 0 : } else Throw( "RefMode not implemented" );
[ - - ][ - - ]
[ - - ][ - - ]
567 : :
568 : : // Communicate extra edges
569 : 177 : comExtra();
570 : 177 : }
571 : :
572 : : void
573 [ + + ]: 181 : Refiner::comExtra()
574 : : // *****************************************************************************
575 : : // Communicate extra edges along chare boundaries
576 : : // *****************************************************************************
577 : : {
578 : : // Export extra added nodes on our mesh chunk boundary to other chares
579 [ + + ]: 181 : if (m_ch.empty()) {
580 : 41 : correctref();
581 : : } else {
582 [ + + ]: 456 : for (auto c : m_ch) { // for all chares we share at least an edge with
583 [ + - ][ - + ]: 948 : thisProxy[c].addRefBndEdges(thisIndex, m_localEdgeData, m_intermediates);
[ - - ]
584 : : }
585 : : }
586 : 181 : }
587 : :
588 : : void
589 : 316 : Refiner::addRefBndEdges(
590 : : int fromch,
591 : : const AMR::EdgeData& ed,
592 : : const std::unordered_set< std::size_t >& intermediates )
593 : : // *****************************************************************************
594 : : //! Receive edges on our chare boundary from other chares
595 : : //! \param[in] fromch Chare call coming from
596 : : //! \param[in] ed Edges on chare boundary
597 : : //! \param[in] intermediates Intermediate nodes
598 : : //! \details Other than update remoteEdge data, this function also updates
599 : : //! locking information for such edges whos nodes are marked as intermediate
600 : : //! by neighboring chares.
601 : : // *****************************************************************************
602 : : {
603 : : // Save/augment buffers of edge data for each sender chare
604 : : auto& red = m_remoteEdgeData[ fromch ];
605 : : auto& re = m_remoteEdges[ fromch ];
606 : : using edge_data_t = std::tuple< Edge, int, int, AMR::Edge_Lock_Case >;
607 [ + + ]: 374792 : for (const auto& [ edge, data ] : ed) {
608 : 374476 : red.push_back( edge_data_t{ edge, std::get<0>(data), std::get<1>(data),
609 : : std::get<2>(data) } );
610 : 374476 : re.push_back( edge );
611 : : }
612 : :
613 : : // Add intermediates to mesh refiner lib
614 : : // needs to be done only when mesh has been actually updated, i.e. first iter
615 [ + + ]: 316 : if (m_ncit == 0) {
616 : 624 : auto esup = tk::genEsup( m_inpoel, 4 );
617 [ + - ]: 624 : auto psup = tk::genPsup( m_inpoel, 4, esup );
618 [ + + ]: 406 : for (const auto g : intermediates) {
619 : 94 : auto l = m_lid.find( g ); // convert to local node ids
620 [ + + ]: 94 : if (l != end(m_lid)) {
621 : : // lock all edges connected to intermediate node
622 : 2 : auto p = l->second;
623 [ + + ]: 30 : for (auto q : tk::Around(psup,p)) {
624 [ + + ]: 28 : AMR::edge_t e(m_rid[p], m_rid[q]);
625 [ + - ]: 28 : auto& refedge = m_refiner.tet_store.edge_store.get(e);
626 [ + + ]: 28 : if (refedge.lock_case == AMR::Edge_Lock_Case::unlocked) {
627 : 1 : refedge.lock_case = AMR::Edge_Lock_Case::temporary; //intermediate;
628 : 1 : refedge.needs_refining = 0;
629 : : }
630 : : }
631 : : }
632 : : }
633 : : }
634 : :
635 : : // Heard from every worker we share at least a single edge with
636 [ + + ]: 316 : if (++m_nref == m_ch.size()) {
637 : 140 : m_nref = 0;
638 : :
639 : 140 : updateBndEdgeData();
640 : :
641 : 140 : std::vector< std::size_t > meshdata{ m_meshid };
642 [ + - ]: 140 : contribute( meshdata, CkReduction::max_ulong,
643 : : m_cbr.get< tag::compatibility >() );
644 : : }
645 : 316 : }
646 : :
647 : : void
648 : 181 : Refiner::correctref()
649 : : // *****************************************************************************
650 : : // Correct extra edges to arrive at conforming mesh across chare boundaries
651 : : //! \details This function is called repeatedly until there is not a a single
652 : : //! edge that needs correction for the whole distributed problem to arrive at
653 : : //! a conforming mesh across chare boundaries during a mesh refinement step.
654 : : // *****************************************************************************
655 : : {
656 : : auto unlocked = AMR::Edge_Lock_Case::unlocked;
657 : :
658 : : // Storage for edge data that need correction to yield a conforming mesh
659 : : AMR::EdgeData extra;
660 : : std::size_t neigh_extra(0);
661 : :
662 : : // Vars for debugging purposes
663 : : std::size_t nlocked(0);
664 : : std::array< std::size_t, 4 > ncorrcase{{0,0,0,0}};
665 : :
666 : : // loop through all edges shared with other chares
667 [ + + ]: 497 : for (const auto& [ neighborchare, edgedata ] : m_remoteEdgeData) {
668 : : for (const auto& [edge,remote_needs_refining,remote_needs_derefining,
669 [ + + ]: 374792 : remote_lock_case] : edgedata) {
670 : : // find local data of remote edge
671 : : auto it = m_localEdgeData.find( edge );
672 [ + + ]: 374476 : if (it != end(m_localEdgeData)) {
673 : : auto& local = it->second;
674 : : auto& local_needs_refining = std::get<0>(local);
675 : : auto& local_needs_derefining = std::get<1>(local);
676 : : auto& local_lock_case = std::get<2>(local);
677 : :
678 : 17154 : auto local_needs_refining_orig = local_needs_refining;
679 : 17154 : auto local_needs_derefining_orig = local_needs_derefining;
680 [ + + ]: 17154 : auto local_lock_case_orig = local_lock_case;
681 : :
682 : : Assert( !(local_lock_case > unlocked && local_needs_refining),
683 : : "Invalid local edge: locked & needs refining" );
684 : : Assert( !(remote_lock_case > unlocked && remote_needs_refining),
685 : : "Invalid remote edge: locked & needs refining" );
686 : : Assert( !(local_needs_derefining == 1 && local_needs_refining > 0),
687 : : "Invalid local edge: needs refining and derefining" );
688 : :
689 : : // The parallel compatibility (par-compat) algorithm
690 : :
691 : : // compute lock from local and remote locks as most restrictive
692 : 17154 : local_lock_case = std::max( local_lock_case, remote_lock_case );
693 : :
694 [ + + ]: 17154 : if (local_lock_case > unlocked) {
695 : 524 : local_needs_refining = 0;
696 : : if (local_needs_refining != local_needs_refining_orig ||
697 : : local_lock_case != local_lock_case_orig)
698 : : ++ncorrcase[0];
699 : : }
700 : :
701 : : // Possible combinations of remote-local ref-deref decisions
702 : : // rows 1, 5, 9: no action needed.
703 : : // rows 4, 7, 8: no action on local-side; comm needed.
704 : : //
705 : : // LOCAL | REMOTE | Result
706 : : // 1 d | d | d
707 : : // 2 d | - | -
708 : : // 3 d | r | r
709 : : // 4 - | d | -
710 : : // 5 - | - | -
711 : : // 6 - | r | r
712 : : // 7 r | d | r
713 : : // 8 r | - | r
714 : : // 9 r | r | r
715 : :
716 : : // Rows 3, 6
717 : : // If remote wants to refine
718 [ + + ]: 17154 : if (remote_needs_refining == 1) {
719 [ + + ]: 5191 : if (local_lock_case == unlocked) {
720 : 5170 : local_needs_refining = 1;
721 : 5170 : local_needs_derefining = false;
722 : : if (local_needs_refining != local_needs_refining_orig ||
723 : : local_needs_derefining != local_needs_derefining_orig)
724 : : ++ncorrcase[1];
725 : : }
726 : : else {
727 : : ++nlocked;
728 : : }
729 : : }
730 : :
731 : : // Row 2
732 : : // If remote neither wants to refine nor derefine
733 [ + + ][ + + ]: 17154 : if (remote_needs_refining == 0 && remote_needs_derefining == false) {
734 : 677 : local_needs_derefining = false;
735 : : if (local_needs_derefining != local_needs_derefining_orig)
736 : : ++ncorrcase[2];
737 : : }
738 : :
739 : : // Row 1: special case
740 : : // If remote wants to deref-ref (either of 8:4, 8:2, 4:2)
741 : : // and local does not want to refine (neither pure ref nor deref-ref)
742 [ - + ][ - - ]: 17154 : if (remote_needs_refining == 2 && local_needs_refining == 0) {
743 [ - - ]: 0 : if (local_lock_case == unlocked) {
744 : 0 : local_needs_refining = 1;
745 : 0 : local_needs_derefining = false;
746 : : if (local_needs_refining != local_needs_refining_orig ||
747 : : local_needs_derefining != local_needs_derefining_orig)
748 : : ++ncorrcase[3];
749 : : }
750 : : else {
751 : : ++nlocked;
752 : : }
753 : : }
754 : :
755 : : // Rows 4, 7, 8
756 : :
757 : : // if the remote sent us data that makes us change our local state,
758 : : // e.g., local{-,0} + remote{r,0} -> local{r,0}, i.e., I changed my
759 : : // state I need to tell the world about it
760 [ + + ]: 17154 : if (local_lock_case != local_lock_case_orig ||
761 [ + - ]: 17124 : local_needs_refining != local_needs_refining_orig ||
762 [ - + ]: 17124 : local_needs_derefining != local_needs_derefining_orig)
763 : : {
764 : 60 : auto l1 = tk::cref_find( m_lid, edge[0] );
765 : 30 : auto l2 = tk::cref_find( m_lid, edge[1] );
766 : : Assert( l1 != l2, "Edge end-points local ids are the same" );
767 [ + + ]: 30 : auto r1 = m_rid[ l1 ];
768 [ + + ]: 30 : auto r2 = m_rid[ l2 ];
769 : : Assert( r1 != r2, "Edge end-points refiner ids are the same" );
770 : : //std::cout << thisIndex << ": " << r1 << ", " << r2 << std::endl;
771 : : //if (m_refiner.tet_store.edge_store.get(AMR::edge_t(r1,r2)).lock_case > local_lock_case) {
772 : : // std::cout << thisIndex << ": edge " << r1 << "-" << r2 <<
773 : : // "; prev=" << local_lock_case_orig <<
774 : : // "; new=" << local_lock_case <<
775 : : // "; amr-lib=" << m_refiner.tet_store.edge_store.get(AMR::edge_t(r1,r2)).lock_case <<
776 : : // " | parcompatiter " << m_ncit << std::endl;
777 : : //}
778 [ + + ][ + - ]: 54 : extra[ {{ std::min(r1,r2), std::max(r1,r2) }} ] =
779 : 30 : { local_needs_refining, local_needs_derefining, local_lock_case };
780 : : }
781 : : // or if the remote data is inconsistent with what I think, e.g.,
782 : : // local{r,0} + remote{-,0} -> local{r,0}, i.e., the remote does not
783 : : // yet agree, so another par-compat iteration will be pursued. but
784 : : // I do not have to locally run ref-compat.
785 [ + + ]: 17124 : else if (local_lock_case != remote_lock_case ||
786 [ + - ]: 17094 : local_needs_refining != remote_needs_refining ||
787 [ - + ]: 17094 : local_needs_derefining != remote_needs_derefining)
788 : : {
789 : 30 : ++neigh_extra;
790 : : }
791 : : }
792 : : }
793 : : }
794 : :
795 : : m_remoteEdgeData.clear();
796 : 181 : m_extra = extra.size();
797 : : //std::cout << thisIndex << ": amr correction reqd for nedge: " << m_extra << std::endl;
798 : : //std::cout << thisIndex << ": amr correction reqd for neighbor edges: " << neigh_extra << std::endl;
799 : : //std::cout << thisIndex << ": edge counts by correction case: " << ncorrcase[0]
800 : : // << ", " << ncorrcase[1] << ", " << ncorrcase[2] << ", " << ncorrcase[3] << std::endl;
801 : : //std::cout << thisIndex << ": locked edges that req corr: " << nlocked << std::endl;
802 : :
803 [ + + ]: 181 : if (!extra.empty()) {
804 : : //std::cout << thisIndex << ": redoing markings" << std::endl;
805 : : // Do refinement compatibility (ref-compat) for edges that need correction
806 [ + - ]: 4 : m_refiner.mark_error_refinement_corr( extra );
807 : 4 : ++m_ncit;
808 : : // Update our extra-edge store based on refiner
809 [ + - ]: 4 : updateEdgeData();
810 : : m_remoteEdges.clear();
811 : : }
812 [ + - ]: 177 : else if (neigh_extra == 0) {
813 : 177 : m_ncit = 0;
814 : : }
815 : :
816 : : // Aggregate number of extra edges that still need correction and some
817 : : // refinement/derefinement statistics
818 : : const auto& tet_store = m_refiner.tet_store;
819 : : std::vector< std::size_t >
820 : 181 : m{ m_meshid,
821 : 181 : m_extra,
822 : : tet_store.marked_refinements.size(),
823 : : tet_store.marked_derefinements.size(),
824 [ + - ]: 181 : static_cast< std::underlying_type_t< RefMode > >( m_mode ) };
825 [ + - ]: 181 : contribute( m, CkReduction::sum_ulong, m_cbr.get< tag::matched >() );
826 : 181 : }
827 : :
828 : : void
829 : 358 : Refiner::updateEdgeData()
830 : : // *****************************************************************************
831 : : // Query AMR lib and update our local store of edge data
832 : : // *****************************************************************************
833 : : {
834 : : m_localEdgeData.clear();
835 : : m_intermediates.clear();
836 : :
837 : : // This currently takes ALL edges from the AMR lib, i.e., on the whole
838 : : // domain. We should eventually only collect edges here that are shared with
839 : : // other chares.
840 : : const auto& ref_edges = m_refiner.tet_store.edge_store.edges;
841 : 358 : const auto& refinpoel = m_refiner.tet_store.get_active_inpoel();
842 : :
843 [ + + ]: 589912 : for (std::size_t e=0; e<refinpoel.size()/4; ++e) {
844 : 589554 : auto A = refinpoel[e*4+0];
845 : 589554 : auto B = refinpoel[e*4+1];
846 : 589554 : auto C = refinpoel[e*4+2];
847 : 589554 : auto D = refinpoel[e*4+3];
848 : : std::array<Edge,6> edges{{ {{A,B}}, {{B,C}}, {{A,C}},
849 : 589554 : {{A,D}}, {{B,D}}, {{C,D}} }};
850 [ + + ]: 4126878 : for (const auto& ed : edges) {
851 [ + + ]: 5772966 : auto ae = AMR::edge_t{{{ std::min(ed[0],ed[1]), std::max(ed[0],ed[1]) }}};
852 : 3537324 : auto r = tk::cref_find( ref_edges, ae );
853 : 7074648 : const auto ged = Edge{{ m_gid[ tk::cref_find( m_lref, ed[0] ) ],
854 : 7074648 : m_gid[ tk::cref_find( m_lref, ed[1] ) ] }};
855 : : m_localEdgeData[ ged ] = { r.needs_refining, r.needs_derefining,
856 : : r.lock_case };
857 : : }
858 : : }
859 : :
860 : : // Build intermediates to send. This currently takes ALL intermediates from
861 : : // the AMR lib, i.e., on the whole domain. We should eventually only collect
862 : : // edges here that are shared with other chares.
863 [ + + ]: 728 : for (const auto& r : m_refiner.tet_store.intermediate_list) {
864 : 740 : m_intermediates.insert( m_gid[ tk::cref_find( m_lref, r ) ] );
865 : : }
866 : 358 : }
867 : :
868 : : void
869 : 140 : Refiner::updateBndEdgeData()
870 : : // *****************************************************************************
871 : : // Query AMR lib and update our local store of boundary edge data
872 : : // *****************************************************************************
873 : : {
874 : : // This currently takes ALL edges from the AMR lib, i.e., on the whole
875 : : // domain. We should eventually only collect edges here that are shared with
876 : : // other chares.
877 : : const auto& ref_edges = m_refiner.tet_store.edge_store.edges;
878 : 140 : const auto& refinpoel = m_refiner.tet_store.get_active_inpoel();
879 : :
880 [ + + ]: 109715 : for (std::size_t e=0; e<refinpoel.size()/4; ++e) {
881 : 109575 : auto A = refinpoel[e*4+0];
882 : 109575 : auto B = refinpoel[e*4+1];
883 : 109575 : auto C = refinpoel[e*4+2];
884 : 109575 : auto D = refinpoel[e*4+3];
885 : : std::array<Edge,6> edges{{ {{A,B}}, {{B,C}}, {{A,C}},
886 : 109575 : {{A,D}}, {{B,D}}, {{C,D}} }};
887 [ + + ]: 767025 : for (const auto& ed : edges) {
888 [ + + ]: 1076269 : auto ae = AMR::edge_t{{{ std::min(ed[0],ed[1]), std::max(ed[0],ed[1]) }}};
889 : 657450 : auto r = tk::cref_find( ref_edges, ae );
890 : 1314900 : const auto ged = Edge{{ m_gid[ tk::cref_find( m_lref, ed[0] ) ],
891 : 1314900 : m_gid[ tk::cref_find( m_lref, ed[1] ) ] }};
892 : : // only update edges that are on chare boundary OR unlocked
893 [ + - ]: 657450 : if (m_localEdgeData.find(ged) == m_localEdgeData.end() ||
894 [ + + ]: 657450 : std::get<2>(m_localEdgeData[ged]) == AMR::Edge_Lock_Case::unlocked) {
895 : : m_localEdgeData[ ged ] = { r.needs_refining, r.needs_derefining,
896 : : r.lock_case };
897 : : }
898 : : }
899 : : }
900 : 140 : }
901 : :
902 : : std::tuple< std::vector< std::string >,
903 : : std::vector< std::vector< tk::real > >,
904 : : std::vector< std::string >,
905 : : std::vector< std::vector< tk::real > > >
906 : 149 : Refiner::refinementFields() const
907 : : // *****************************************************************************
908 : : // Collect mesh output fields from refiner lib
909 : : //! \return Names and fields of mesh refinement data in mesh cells and nodes
910 : : // *****************************************************************************
911 : : {
912 : : // Find number of nodes in current mesh
913 : 149 : auto npoin = tk::npoin_in_graph( m_inpoel );
914 : : // Generate edges surrounding points in current mesh
915 : 298 : auto esup = tk::genEsup( m_inpoel, 4 );
916 : :
917 : : // Update solution on current mesh
918 [ + - ]: 149 : const auto& u = solution( npoin, esup );
919 : : Assert( u.nunk() == npoin, "Solution uninitialized or wrong size" );
920 : :
921 : : // Compute error in edges on current mesh
922 [ + - ]: 149 : auto edgeError = errorsInEdges( npoin, esup, u );
923 : :
924 : : // Transfer error from edges to cells for field output
925 [ + - ]: 149 : std::vector< tk::real > error( m_inpoel.size()/4, 0.0 );
926 [ + + ]: 226755 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
927 : 226606 : auto A = m_inpoel[e*4+0];
928 : 226606 : auto B = m_inpoel[e*4+1];
929 : 226606 : auto C = m_inpoel[e*4+2];
930 : 226606 : auto D = m_inpoel[e*4+3];
931 : : std::array<Edge,6> edges{{ {{A,B}}, {{B,C}}, {{A,C}},
932 : 226606 : {{A,D}}, {{B,D}}, {{C,D}} }};
933 : : // sum error from edges to elements
934 [ + + ]: 1586242 : for (const auto& ed : edges) error[e] += tk::cref_find( edgeError, ed );
935 : 226606 : error[e] /= 6.0; // assign edge-average error to element
936 : : }
937 : :
938 : : // Prepare element fields with mesh refinement data
939 : : std::vector< std::string >
940 [ + - ][ + - ]: 298 : elemfieldnames{ "refinement level", "cell type", "error" };
[ + - ][ - + ]
[ - - ]
941 : 149 : auto& tet_store = m_refiner.tet_store;
942 : : std::vector< std::vector< tk::real > > elemfields{
943 : : tet_store.get_refinement_level_list(),
944 : : tet_store.get_cell_type_list(),
945 [ + - ][ + - ]: 298 : error };
[ + - ][ - + ]
[ + - ][ - - ]
946 : :
947 : : using tuple_t = std::tuple< std::vector< std::string >,
948 : : std::vector< std::vector< tk::real > >,
949 : : std::vector< std::string >,
950 : : std::vector< std::vector< tk::real > > >;
951 : 298 : return tuple_t{ elemfieldnames, elemfields, {}, {} };
952 : : }
953 : :
954 : : void
955 : 149 : Refiner::writeMesh( const std::string& basefilename,
956 : : uint64_t itr,
957 : : tk::real t,
958 : : CkCallback c ) const
959 : : // *****************************************************************************
960 : : // Output mesh to file(s)
961 : : //! \param[in] basefilename File name to append to
962 : : //! \param[in] itr Iteration count since a new mesh
963 : : //! \param[in] t "Physical time" to write to file. "Time" here is used to
964 : : //! designate a new time step at which the mesh is saved.
965 : : //! \param[in] c Function to continue with after the write
966 : : // *****************************************************************************
967 : : {
968 : 149 : auto r = refinementFields();
969 : : auto& elemfieldnames = std::get< 0 >( r );
970 : : auto& elemfields = std::get< 1 >( r );
971 : : auto& nodefieldnames = std::get< 2 >( r );
972 : : auto& nodefields = std::get< 3 >( r );
973 : :
974 : : // Prepare solution field names: depvar + component id for all eqs
975 : 149 : auto nprop = g_inputdeck.get< tag::ncomp >();
976 : 149 : std::vector< std::string > solfieldnames;
977 [ + + ]: 346 : for (std::size_t i=0; i<nprop; ++i) {
978 : : solfieldnames.push_back(
979 [ + - ][ + - ]: 394 : g_inputdeck.get< tag::depvar >()[0] + std::to_string(i+1) );
[ - + ][ - - ]
980 : : }
981 : : Assert( solfieldnames.size() == nprop, "Size mismatch" );
982 : :
983 : 149 : const auto scheme = g_inputdeck.get< tag::scheme >();
984 [ + - ][ + - ]: 149 : const auto centering = ctr::Scheme().centering( scheme );
985 [ + + ]: 149 : auto t0 = g_inputdeck.get< tag::t0 >();
986 : :
987 : : // list of nodes/elements at which box ICs are defined
988 : 149 : std::vector< std::unordered_set< std::size_t > > inbox;
989 : : tk::real V = 1.0;
990 : : std::vector< tk::real > blkvols;
991 : : std::unordered_map< std::size_t, std::set< std::size_t > > nodeblockid,
992 : : elemblockid;
993 : :
994 : : // Prepare node or element fields for output to file
995 [ + + ]: 149 : if (centering == tk::Centering::NODE) {
996 : :
997 : : // Augment element field names with solution variable names + field ids
998 : : nodefieldnames.insert( end(nodefieldnames),
999 [ + - ]: 83 : begin(solfieldnames), end(solfieldnames) );
1000 : :
1001 : : // Evaluate initial conditions on current mesh at t0
1002 [ + - ]: 83 : tk::Fields u( m_coord[0].size(), nprop );
1003 [ + - ]: 83 : g_cgpde[m_meshid].initialize( m_coord, u, t0, V, inbox, blkvols,
1004 : : nodeblockid );
1005 : :
1006 : : // Extract all scalar components from solution for output to file
1007 [ + + ]: 198 : for (std::size_t i=0; i<nprop; ++i)
1008 [ + - ][ - - ]: 230 : nodefields.push_back( u.extract_comp( i ) );
1009 : :
1010 [ + - ]: 66 : } else if (centering == tk::Centering::ELEM) {
1011 : :
1012 : : // Augment element field names with solution variable names + field ids
1013 : : elemfieldnames.insert( end(elemfieldnames),
1014 [ + - ]: 66 : begin(solfieldnames), end(solfieldnames) );
1015 : :
1016 : 66 : auto ndof = g_inputdeck.get< tag::ndof >();
1017 [ + - ]: 66 : tk::Fields lhs( m_inpoel.size()/4, ndof*nprop );
1018 : :
1019 : : // Generate left hand side for DG and evaluate initial conditions on
1020 : : // current mesh at t0
1021 [ + - ]: 66 : auto geoElem = tk::genGeoElemTet( m_inpoel, m_coord );
1022 : : auto u = lhs;
1023 [ - + ]: 66 : if (scheme == ctr::SchemeType::FV) {
1024 [ - - ]: 0 : g_fvpde[m_meshid].lhs( geoElem, lhs );
1025 [ - - ]: 0 : g_fvpde[m_meshid].initialize( lhs, m_inpoel, m_coord, inbox, elemblockid,
1026 [ - - ]: 0 : u, t0, m_inpoel.size()/4 );
1027 : : }
1028 : : else {
1029 [ + - ]: 66 : g_dgpde[m_meshid].lhs( geoElem, lhs );
1030 [ + - ]: 66 : g_dgpde[m_meshid].initialize( lhs, m_inpoel, m_coord, inbox, elemblockid,
1031 [ + - ]: 66 : u, t0, m_inpoel.size()/4 );
1032 : : }
1033 : :
1034 : : // Extract all scalar components from solution for output to file
1035 [ + + ]: 148 : for (std::size_t i=0; i<nprop; ++i)
1036 [ + - ][ - - ]: 164 : elemfields.push_back( u.extract_comp( i ) );
1037 : : }
1038 : :
1039 : : // Output mesh
1040 : : m_meshwriter[ CkNodeFirst( CkMyNode() ) ].
1041 [ + - ][ - + ]: 447 : write( m_meshid, /*meshoutput = */ true, /*fieldoutput = */ true, itr, 1, t,
[ - - ]
1042 [ + - ]: 149 : thisIndex, basefilename, m_inpoel, m_coord, m_bface,
1043 [ + - ][ + - ]: 372 : tk::remap(m_bnode,m_lid), tk::remap(m_triinpoel,m_lid),
[ + + ][ - - ]
1044 : : elemfieldnames, nodefieldnames, {}, {}, elemfields, nodefields, {},
1045 : : {}, {}, c );
1046 : 149 : }
1047 : :
1048 : : void
1049 : 177 : Refiner::perform()
1050 : : // *****************************************************************************
1051 : : // Perform mesh refinement and decide how to continue
1052 : : //! \details First the mesh refiner object is called to perform a single step
1053 : : //! of mesh refinement. Then, if this function is called during a step
1054 : : //! (potentially multiple levels of) initial AMR, it evaluates whether to do
1055 : : //! another one. If it is called during time stepping, this concludes the
1056 : : //! single mesh refinement step and the new mesh is sent to the PDE worker
1057 : : //! (Discretization).
1058 : : // *****************************************************************************
1059 : : {
1060 : : // Save old tets and their ids before performing refinement. Outref is always
1061 : : // followed by outderef, so to the outside world, the mesh is uchanged, thus
1062 : : // no update.
1063 [ + + ]: 177 : if (m_mode != RefMode::OUTREF && m_mode != RefMode::OUTDEREF) {
1064 : : m_oldTets.clear();
1065 [ + + ]: 129885 : for (const auto& [ id, tet ] : m_refiner.tet_store.tets) {
1066 : : m_oldTets.insert( tet );
1067 : : }
1068 : 111 : m_oldntets = m_oldTets.size();
1069 : : }
1070 : :
1071 [ + + ]: 177 : if (m_mode == RefMode::T0REF) {
1072 : :
1073 : : // Refine mesh based on next initial refinement type
1074 [ + - ]: 111 : if (!m_initref.empty()) {
1075 : 111 : auto r = m_initref.back(); // consume (reversed) list from its back
1076 [ + + ]: 111 : if (r == ctr::AMRInitialType::UNIFORM_DEREFINE)
1077 : 36 : m_refiner.perform_derefinement();
1078 : : else
1079 : 75 : m_refiner.perform_refinement();
1080 : : }
1081 : :
1082 : : } else {
1083 : :
1084 : : // TODO: does not work yet, fix as above
1085 : 66 : m_refiner.perform_refinement();
1086 : 66 : m_refiner.perform_derefinement();
1087 : : }
1088 : :
1089 : : // Remove temporary edge-locks resulting from the parallel compatibility
1090 : 177 : m_refiner.remove_edge_locks(1);
1091 : 177 : m_refiner.remove_edge_temp_locks();
1092 : :
1093 : : //auto& tet_store = m_refiner.tet_store;
1094 : : //std::cout << "before ref: " << tet_store.marked_refinements.size() << ", " << tet_store.marked_derefinements.size() << ", " << tet_store.size() << ", " << tet_store.get_active_inpoel().size() << '\n';
1095 : : //std::cout << "after ref: " << tet_store.marked_refinements.size() << ", " << tet_store.marked_derefinements.size() << ", " << tet_store.size() << ", " << tet_store.get_active_inpoel().size() << '\n';
1096 : : //std::cout << "after deref: " << tet_store.marked_refinements.size() << ", " << tet_store.marked_derefinements.size() << ", " << tet_store.size() << ", " << tet_store.get_active_inpoel().size() << '\n';
1097 : :
1098 : : // Update volume and boundary mesh
1099 : 177 : updateMesh();
1100 : :
1101 : : // Save mesh at every initial refinement step (mainly for debugging). Will
1102 : : // replace with just a 'next()' in production.
1103 [ + + ]: 177 : if (m_mode == RefMode::T0REF) {
1104 : :
1105 : 111 : auto l = m_ninitref - m_initref.size() + 1; // num initref steps completed
1106 : 111 : auto t0 = g_inputdeck.get< tag::t0 >();
1107 : : // Generate times equally subdividing t0-1...t0 to ninitref steps
1108 : 111 : auto t =
1109 : 111 : t0 - 1.0 + static_cast<tk::real>(l)/static_cast<tk::real>(m_ninitref);
1110 : : auto itr = l;
1111 : : // Output mesh after refinement step
1112 [ + - ][ + - ]: 222 : writeMesh( "t0ref", itr, t,
1113 [ + - ][ - + ]: 333 : CkCallback( CkIndex_Refiner::next(), thisProxy[thisIndex] ) );
[ - - ]
1114 : :
1115 : : } else {
1116 : :
1117 : 66 : next();
1118 : :
1119 : : }
1120 : 177 : }
1121 : :
1122 : : void
1123 : 177 : Refiner::next()
1124 : : // *****************************************************************************
1125 : : // Continue after finishing a refinement step
1126 : : // *****************************************************************************
1127 : : {
1128 [ + + ]: 177 : if (m_mode == RefMode::T0REF) {
1129 : :
1130 : : // Remove initial mesh refinement step from list
1131 [ + - ]: 111 : if (!m_initref.empty()) m_initref.pop_back();
1132 : : // Continue to next initial AMR step or finish
1133 [ + + ]: 111 : if (!m_initref.empty()) t0ref(); else endt0ref();
1134 : :
1135 [ - + ]: 66 : } else if (m_mode == RefMode::DTREF) {
1136 : :
1137 : : // Send new mesh, solution, and communication data back to PDE worker
1138 : 0 : m_scheme[m_meshid].ckLocal< Scheme::resizePostAMR >( thisIndex, m_ginpoel,
1139 : 0 : m_el, m_coord, m_addedNodes, m_addedTets, m_removedNodes, m_amrNodeMap,
1140 : 0 : m_nodeCommMap, m_bface, m_bnode, m_triinpoel, m_elemblockid );
1141 : :
1142 [ + + ]: 66 : } else if (m_mode == RefMode::OUTREF) {
1143 : :
1144 : : // Store field output mesh
1145 : 33 : m_outref_ginpoel = m_ginpoel;
1146 : : m_outref_el = m_el;
1147 : : m_outref_coord = m_coord;
1148 : : m_outref_addedNodes = m_addedNodes;
1149 : : m_outref_addedTets = m_addedTets;
1150 : : m_outref_nodeCommMap = m_nodeCommMap;
1151 : : m_outref_bface = m_bface;
1152 : : m_outref_bnode = m_bnode;
1153 : 33 : m_outref_triinpoel = m_triinpoel;
1154 : :
1155 : : // Derefine mesh to the state previous to field output
1156 [ + - ]: 66 : outref( m_outref_bface, m_outref_bnode, m_outref_triinpoel, m_writeCallback,
1157 : : RefMode::OUTDEREF );
1158 : :
1159 [ + - ]: 33 : } else if (m_mode == RefMode::OUTDEREF) {
1160 : :
1161 : : // Send field output mesh to PDE worker
1162 : 66 : m_scheme[m_meshid].ckLocal< Scheme::extractFieldOutput >( thisIndex,
1163 : 33 : m_outref_ginpoel, m_outref_el, m_outref_coord, m_outref_addedNodes,
1164 : 33 : m_outref_addedTets, m_outref_nodeCommMap, m_outref_bface, m_outref_bnode,
1165 : 33 : m_outref_triinpoel, m_writeCallback );
1166 : :
1167 [ - - ][ - - ]: 0 : } else Throw( "RefMode not implemented" );
[ - - ][ - - ]
[ - - ][ - - ]
1168 : 177 : }
1169 : :
1170 : : void
1171 : 1943 : Refiner::endt0ref()
1172 : : // *****************************************************************************
1173 : : // Finish initial mesh refinement
1174 : : //! \details This function is called as after initial mesh refinement has
1175 : : //! finished. If initial mesh reifnement was not configured by the user, this
1176 : : //! is the point where we continue after the constructor, by computing the
1177 : : //! total number of elements across the whole problem.
1178 : : // *****************************************************************************
1179 : : {
1180 : : // create sorter Charm++ chare array elements using dynamic insertion
1181 [ + - ]: 3886 : m_sorter[ thisIndex ].insert( m_meshid, m_host, m_meshwriter, m_cbs, m_scheme,
1182 [ + - ][ - + ]: 3886 : CkCallback(CkIndex_Refiner::reorder(), thisProxy[thisIndex]), m_ginpoel,
[ - + ][ - - ]
[ - - ]
1183 [ + - ]: 1943 : m_coordmap, m_el, m_bface, m_triinpoel, m_bnode, m_elemblockid, m_nchare );
1184 : :
1185 : : // Compute final number of cells across whole problem
1186 : : std::vector< std::size_t >
1187 : 1943 : meshdata{ m_meshid, m_ginpoel.size()/4, m_coord[0].size() };
1188 [ + - ]: 1943 : contribute( meshdata, CkReduction::sum_ulong, m_cbr.get< tag::refined >() );
1189 : :
1190 : : // // Free up memory if no dtref
1191 : : // if (!g_inputdeck.get< tag::amr, tag::dtref >()) {
1192 : : // tk::destroy( m_ginpoel );
1193 : : // tk::destroy( m_el );
1194 : : // tk::destroy( m_coordmap );
1195 : : // tk::destroy( m_coord );
1196 : : // tk::destroy( m_bface );
1197 : : // tk::destroy( m_bnode );
1198 : : // tk::destroy( m_triinpoel );
1199 : : // tk::destroy( m_initref );
1200 : : // tk::destroy( m_ch );
1201 : : // tk::destroy( m_edgech );
1202 : : // tk::destroy( m_chedge );
1203 : : // tk::destroy( m_localEdgeData );
1204 : : // tk::destroy( m_remoteEdgeData );
1205 : : // tk::destroy( m_remoteEdges );
1206 : : // tk::destroy( m_intermediates );
1207 : : // tk::destroy( m_nodeCommMap );
1208 : : // tk::destroy( m_oldTets );
1209 : : // tk::destroy( m_addedNodes );
1210 : : // tk::destroy( m_addedTets );
1211 : : // tk::destroy( m_coarseBndFaces );
1212 : : // tk::destroy( m_coarseBndNodes );
1213 : : // tk::destroy( m_rid );
1214 : : // // tk::destroy( m_oldrid );
1215 : : // tk::destroy( m_lref );
1216 : : // // tk::destroy( m_oldparent );
1217 : : // }
1218 : 1943 : }
1219 : :
1220 : : void
1221 : 90 : Refiner::uniformRefine()
1222 : : // *****************************************************************************
1223 : : // Do uniform mesh refinement
1224 : : // *****************************************************************************
1225 : : {
1226 : : // Do uniform refinement
1227 : 90 : m_refiner.mark_uniform_refinement();
1228 : :
1229 : : // Update our extra-edge store based on refiner
1230 : 90 : updateEdgeData();
1231 : :
1232 : : // Set number of extra edges to be zero, skipping correction (if this is the
1233 : : // only step in this refinement step)
1234 : 90 : m_extra = 0;
1235 : 90 : }
1236 : :
1237 : : void
1238 : 69 : Refiner::uniformDeRefine()
1239 : : // *****************************************************************************
1240 : : // Do uniform mesh derefinement
1241 : : // *****************************************************************************
1242 : : {
1243 : : // Do uniform derefinement
1244 : 69 : m_refiner.mark_uniform_derefinement();
1245 : :
1246 : : // Update our extra-edge store based on refiner
1247 : 69 : updateEdgeData();
1248 : :
1249 : : // Set number of extra edges to be zero, skipping correction (if this is the
1250 : : // only step in this refinement step)
1251 : 69 : m_extra = 0;
1252 : 69 : }
1253 : :
1254 : : Refiner::EdgeError
1255 : 165 : Refiner::errorsInEdges(
1256 : : std::size_t npoin,
1257 : : const std::pair< std::vector<std::size_t>, std::vector<std::size_t> >& esup,
1258 : : const tk::Fields& u ) const
1259 : : // *****************************************************************************
1260 : : // Compute errors in edges
1261 : : //! \param[in] npoin Number nodes in current mesh (partition)
1262 : : //! \param[in] esup Elements surrounding points linked vectors
1263 : : //! \param[in] u Solution evaluated at mesh nodes for all scalar components
1264 : : //! \return A map associating errors (real values between 0.0 and 1.0 incusive)
1265 : : //! to edges (2 local node IDs)
1266 : : // *****************************************************************************
1267 : : {
1268 : : // Get the indices (in the system of systems) of refinement variables and the
1269 : : // error indicator configured
1270 : 165 : auto ncomp = g_inputdeck.get< tag::ncomp >();
1271 : 165 : auto errtype = g_inputdeck.get< tag::amr, tag::error >();
1272 : :
1273 : : // Compute points surrounding points
1274 : 330 : auto psup = tk::genPsup( m_inpoel, 4, esup );
1275 : :
1276 : : // Compute errors in ICs and define refinement criteria for edges
1277 : : AMR::Error error;
1278 : : EdgeError edgeError;
1279 : :
1280 [ + + ]: 60985 : for (std::size_t p=0; p<npoin; ++p) { // for all mesh nodes on this chare
1281 [ + + ][ + + ]: 721982 : for (auto q : tk::Around(psup,p)) { // for all nodes surrounding p
1282 : : tk::real cmax = 0.0;
1283 : 661162 : AMR::edge_t e(p,q);
1284 [ + + ]: 1550636 : for (std::size_t i=0; i<ncomp; ++i) { // for all refinement variables
1285 [ + - ]: 889474 : auto c = error.scalar( u, e, i, m_coord, m_inpoel, esup, errtype );
1286 [ + + ]: 889474 : if (c > cmax) cmax = c; // find max error at edge
1287 : : }
1288 [ + - ]: 661162 : edgeError[ {{p,q}} ] = cmax; // associate error to edge
1289 : : }
1290 : : }
1291 : :
1292 : 165 : return edgeError;
1293 : : }
1294 : :
1295 : : tk::Fields
1296 [ + - ]: 165 : Refiner::solution( std::size_t npoin,
1297 : : const std::pair< std::vector< std::size_t >,
1298 : : std::vector< std::size_t > >& esup ) const
1299 : : // *****************************************************************************
1300 : : // Update (or evaluate) solution on current mesh
1301 : : //! \param[in] npoin Number nodes in current mesh (partition)
1302 : : //! \param[in] esup Elements surrounding points linked vectors
1303 : : //! \return Solution updated/evaluated for all scalar components
1304 : : // *****************************************************************************
1305 : : {
1306 : : // Get solution whose error to evaluate
1307 : : tk::Fields u;
1308 : :
1309 [ + - ]: 165 : if (m_mode == RefMode::T0REF) {
1310 : :
1311 : : // Evaluate initial conditions at mesh nodes
1312 [ + - ]: 330 : u = nodeinit( npoin, esup );
1313 : :
1314 [ - - ]: 0 : } else if (m_mode == RefMode::DTREF) {
1315 : :
1316 : : // Query current solution
1317 [ - - ]: 0 : u = m_scheme[m_meshid].ckLocal< Scheme::solution >( thisIndex );
1318 : :
1319 : 0 : const auto scheme = g_inputdeck.get< tag::scheme >();
1320 [ - - ][ - - ]: 0 : const auto centering = ctr::Scheme().centering( scheme );
[ - - ]
1321 [ - - ]: 0 : if (centering == tk::Centering::ELEM) {
1322 : :
1323 : : // ...
1324 [ - - ][ - - ]: 0 : Throw("Element-based solution handling not implemented in Refiner");
[ - - ][ - - ]
[ - - ][ - - ]
1325 : :
1326 : : }
1327 : :
1328 [ - - ]: 0 : } else if (m_mode == RefMode::OUTREF) {
1329 : :
1330 : :
1331 : :
1332 [ - - ]: 0 : } else if (m_mode == RefMode::OUTDEREF) {
1333 : :
1334 : :
1335 : :
1336 [ - - ][ - - ]: 0 : } else Throw( "RefMode not implemented" );
[ - - ][ - - ]
[ - - ][ - - ]
1337 : :
1338 : 165 : return u;
1339 : : }
1340 : :
1341 : : void
1342 : 16 : Refiner::errorRefine()
1343 : : // *****************************************************************************
1344 : : // Do error-based mesh refinement and derefinement
1345 : : // *****************************************************************************
1346 : : {
1347 : : // Find number of nodes in old mesh
1348 : 16 : auto npoin = tk::npoin_in_graph( m_inpoel );
1349 : : // Generate edges surrounding points in old mesh
1350 : 32 : auto esup = tk::genEsup( m_inpoel, 4 );
1351 : :
1352 : : // Update solution on current mesh
1353 [ + - ]: 16 : const auto& u = solution( npoin, esup );
1354 : : Assert( u.nunk() == npoin, "Solution uninitialized or wrong size" );
1355 : :
1356 : : using AMR::edge_t;
1357 : : using AMR::edge_tag;
1358 : :
1359 : : // Compute error in edges. Tag edge for refinement if error exceeds
1360 : : // refinement tolerance, tag edge for derefinement if error is below
1361 : : // derefinement tolerance.
1362 : 16 : auto tolref = g_inputdeck.get< tag::amr, tag::tol_refine >();
1363 [ + - ]: 16 : auto tolderef = g_inputdeck.get< tag::amr, tag::tol_derefine >();
1364 : : std::vector< std::pair< edge_t, edge_tag > > tagged_edges;
1365 [ + - ][ + + ]: 3394 : for (const auto& e : errorsInEdges(npoin,esup,u)) {
1366 [ + + ]: 3362 : if (e.second > tolref) {
1367 [ + + ][ + - ]: 3312 : tagged_edges.push_back( { edge_t( m_rid[e.first[0]], m_rid[e.first[1]] ),
1368 : : edge_tag::REFINE } );
1369 [ + + ]: 1706 : } else if (e.second < tolderef) {
1370 [ + + ][ + - ]: 3132 : tagged_edges.push_back( { edge_t( m_rid[e.first[0]], m_rid[e.first[1]] ),
1371 : : edge_tag::DEREFINE } );
1372 : : }
1373 : : }
1374 : :
1375 : : // Do error-based refinement
1376 [ + - ]: 16 : m_refiner.mark_error_refinement( tagged_edges );
1377 : :
1378 : : // Update our extra-edge store based on refiner
1379 [ + - ]: 16 : updateEdgeData();
1380 : :
1381 : : // Set number of extra edges to a nonzero number, triggering correction
1382 [ + - ]: 16 : m_extra = 1;
1383 : 16 : }
1384 : :
1385 : : void
1386 [ - - ]: 0 : Refiner::edgelistRefine()
1387 : : // *****************************************************************************
1388 : : // Do mesh refinement based on user explicitly tagging edges
1389 : : // *****************************************************************************
1390 : : {
1391 : : // Get user-defined node-pairs (edges) to tag for refinement
1392 : : const auto& edgenodelist = g_inputdeck.get< tag::amr, tag::edgelist >();
1393 : :
1394 [ - - ]: 0 : if (!edgenodelist.empty()) { // if user explicitly tagged edges
1395 : : // Find number of nodes in old mesh
1396 : 0 : auto npoin = tk::npoin_in_graph( m_inpoel );
1397 : : // Generate edges surrounding points in old mesh
1398 : 0 : auto esup = tk::genEsup( m_inpoel, 4 );
1399 [ - - ]: 0 : auto psup = tk::genPsup( m_inpoel, 4, esup );
1400 : :
1401 : : EdgeSet useredges;
1402 [ - - ]: 0 : for (std::size_t i=0; i<edgenodelist.size()/2; ++i)
1403 [ - - ]: 0 : useredges.insert( {{ {edgenodelist[i*2+0], edgenodelist[i*2+1]} }} );
1404 : :
1405 : : using AMR::edge_t;
1406 : : using AMR::edge_tag;
1407 : :
1408 : : // Tag edges the user configured
1409 : : std::vector< std::pair< edge_t, edge_tag > > tagged_edges;
1410 [ - - ]: 0 : for (std::size_t p=0; p<npoin; ++p) // for all mesh nodes on this chare
1411 [ - - ]: 0 : for (auto q : tk::Around(psup,p)) { // for all nodes surrounding p
1412 : 0 : Edge e{{ m_gid[p], m_gid[q] }};
1413 : : auto f = useredges.find(e);
1414 [ - - ]: 0 : if (f != end(useredges)) { // tag edge if on user's list
1415 [ - - ][ - - ]: 0 : tagged_edges.push_back( { edge_t( m_rid[p], m_rid[q] ),
[ - - ]
1416 : : edge_tag::REFINE } );
1417 : : useredges.erase( f );
1418 : : }
1419 : : }
1420 : :
1421 : : // Do error-based refinement
1422 [ - - ]: 0 : m_refiner.mark_error_refinement( tagged_edges );
1423 : :
1424 : : // Update our extra-edge store based on refiner
1425 [ - - ]: 0 : updateEdgeData();
1426 : :
1427 : : // Set number of extra edges to a nonzero number, triggering correction
1428 [ - - ]: 0 : m_extra = 1;
1429 : : }
1430 : 0 : }
1431 : :
1432 : : void
1433 : 2 : Refiner::coordRefine()
1434 : : // *****************************************************************************
1435 : : // Do mesh refinement based on tagging edges based on end-point coordinates
1436 : : // *****************************************************************************
1437 : : {
1438 : : // Get user-defined half-world coordinates
1439 : : const auto& amr_coord = g_inputdeck.get< tag::amr, tag::coords >();
1440 : 2 : auto xminus = amr_coord.get< tag::xminus >();
1441 : 2 : auto xplus = amr_coord.get< tag::xplus >();
1442 : 2 : auto yminus = amr_coord.get< tag::yminus >();
1443 : 2 : auto yplus = amr_coord.get< tag::yplus >();
1444 : 2 : auto zminus = amr_coord.get< tag::zminus >();
1445 : 2 : auto zplus = amr_coord.get< tag::zplus >();
1446 : :
1447 : : // The default is the largest representable double
1448 : : auto eps =
1449 : : std::numeric_limits< tk::real >::epsilon();
1450 : : const auto& amr_defcoord = g_inputdeck_defaults.get< tag::amr, tag::coords >();
1451 : 2 : auto xminus_default = amr_defcoord.get< tag::xminus >();
1452 : 2 : auto xplus_default = amr_defcoord.get< tag::xplus >();
1453 : 2 : auto yminus_default = amr_defcoord.get< tag::yminus >();
1454 : 2 : auto yplus_default = amr_defcoord.get< tag::yplus >();
1455 : 2 : auto zminus_default = amr_defcoord.get< tag::zminus >();
1456 : 2 : auto zplus_default = amr_defcoord.get< tag::zplus >();
1457 : :
1458 : : // Decide if user has configured the half-world
1459 [ - + ]: 2 : bool xm = std::abs(xminus - xminus_default) > eps ? true : false;
1460 : 2 : bool xp = std::abs(xplus - xplus_default) > eps ? true : false;
1461 : 2 : bool ym = std::abs(yminus - yminus_default) > eps ? true : false;
1462 : 2 : bool yp = std::abs(yplus - yplus_default) > eps ? true : false;
1463 : 2 : bool zm = std::abs(zminus - zminus_default) > eps ? true : false;
1464 : 2 : bool zp = std::abs(zplus - zplus_default) > eps ? true : false;
1465 : :
1466 : : using AMR::edge_t;
1467 : : using AMR::edge_tag;
1468 : :
1469 [ - + ][ - - ]: 2 : if (xm || xp || ym || yp || zm || zp) { // if any half-world configured
[ - - ]
1470 : : // Find number of nodes in old mesh
1471 : 2 : auto npoin = tk::npoin_in_graph( m_inpoel );
1472 : : // Generate edges surrounding points in old mesh
1473 : 4 : auto esup = tk::genEsup( m_inpoel, 4 );
1474 [ + - ]: 4 : auto psup = tk::genPsup( m_inpoel, 4, esup );
1475 : : // Get access to node coordinates
1476 : : const auto& x = m_coord[0];
1477 : : const auto& y = m_coord[1];
1478 : : const auto& z = m_coord[2];
1479 : : // Compute edges to be tagged for refinement
1480 : : std::vector< std::pair< edge_t, edge_tag > > tagged_edges;
1481 [ + + ]: 621 : for (std::size_t p=0; p<npoin; ++p) // for all mesh nodes on this chare
1482 [ + + ]: 7181 : for (auto q : tk::Around(psup,p)) { // for all nodes surrounding p
1483 : : Edge e{{p,q}};
1484 : :
1485 : : bool t = true;
1486 [ + - ][ + + ]: 6562 : if (xm) { if (x[p]>xminus && x[q]>xminus) t = false; }
[ + + ]
1487 [ - + ][ - - ]: 6562 : if (xp) { if (x[p]<xplus && x[q]<xplus) t = false; }
[ - - ]
1488 [ - + ][ - - ]: 6562 : if (ym) { if (y[p]>yminus && y[q]>yminus) t = false; }
[ - - ]
1489 [ - + ][ - - ]: 6562 : if (yp) { if (y[p]<yplus && y[q]<yplus) t = false; }
[ - - ]
1490 [ - + ][ - - ]: 6562 : if (zm) { if (z[p]>zminus && z[q]>zminus) t = false; }
[ - - ]
1491 [ - + ][ - - ]: 6562 : if (zp) { if (z[p]<zplus && z[q]<zplus) t = false; }
[ - - ]
1492 : :
1493 [ + + ]: 6562 : if (t) {
1494 [ + + ][ + - ]: 9596 : tagged_edges.push_back( { edge_t( m_rid[e[0]], m_rid[e[1]] ),
[ - - ]
1495 : : edge_tag::REFINE } );
1496 : : }
1497 : : }
1498 : :
1499 : : // Do error-based refinement
1500 [ + - ]: 2 : m_refiner.mark_error_refinement( tagged_edges );
1501 : :
1502 : : // Update our extra-edge store based on refiner
1503 [ + - ]: 2 : updateEdgeData();
1504 : :
1505 : : // Set number of extra edges to a nonzero number, triggering correction
1506 [ + - ]: 2 : m_extra = 1;
1507 : : }
1508 : 2 : }
1509 : :
1510 : : tk::Fields
1511 : 165 : Refiner::nodeinit( std::size_t npoin,
1512 : : const std::pair< std::vector< std::size_t >,
1513 : : std::vector< std::size_t > >& esup ) const
1514 : : // *****************************************************************************
1515 : : // Evaluate initial conditions (IC) at mesh nodes
1516 : : //! \param[in] npoin Number points in mesh (on this chare)
1517 : : //! \param[in] esup Elements surrounding points as linked lists, see tk::genEsup
1518 : : //! \return Initial conditions (evaluated at t0) at nodes
1519 : : // *****************************************************************************
1520 : : {
1521 : 165 : auto t0 = g_inputdeck.get< tag::t0 >();
1522 : 165 : auto nprop = g_inputdeck.get< tag::ncomp >();
1523 : :
1524 : : // Will store nodal ICs
1525 : 165 : tk::Fields u( m_coord[0].size(), nprop );
1526 : :
1527 : : // Evaluate ICs differently depending on nodal or cell-centered discretization
1528 : 165 : const auto scheme = g_inputdeck.get< tag::scheme >();
1529 [ + - ][ + - ]: 330 : const auto centering = ctr::Scheme().centering( scheme );
[ + + ][ - - ]
1530 : : // list of nodes/elements at which box ICs are defined
1531 : 165 : std::vector< std::unordered_set< std::size_t > > inbox;
1532 : : tk::real V = 1.0;
1533 : : std::vector< tk::real > blkvols;
1534 : : std::unordered_map< std::size_t, std::set< std::size_t > > nodeblockid,
1535 : : elemblockid;
1536 : :
1537 [ + + ]: 165 : if (centering == tk::Centering::NODE) {
1538 : :
1539 : : // Evaluate ICs for all scalar components integrated
1540 [ + - ]: 99 : g_cgpde[m_meshid].initialize( m_coord, u, t0, V, inbox, blkvols,
1541 : : nodeblockid );
1542 : :
1543 [ + - ]: 66 : } else if (centering == tk::Centering::ELEM) {
1544 : :
1545 [ + - ]: 66 : auto esuel = tk::genEsuelTet( m_inpoel, esup ); // elems surrounding elements
1546 : : // Initialize cell-based unknowns
1547 [ + - ]: 66 : tk::Fields ue( m_inpoel.size()/4, nprop );
1548 : : auto lhs = ue;
1549 [ + - ]: 66 : auto geoElem = tk::genGeoElemTet( m_inpoel, m_coord );
1550 [ - + ]: 66 : if (scheme == ctr::SchemeType::FV) {
1551 [ - - ]: 0 : g_fvpde[m_meshid].lhs( geoElem, lhs );
1552 [ - - ]: 0 : g_fvpde[m_meshid].initialize( lhs, m_inpoel, m_coord, inbox, elemblockid,
1553 [ - - ]: 0 : ue, t0, esuel.size()/4 );
1554 : : }
1555 : : else {
1556 [ + - ]: 66 : g_dgpde[m_meshid].lhs( geoElem, lhs );
1557 [ + - ]: 66 : g_dgpde[m_meshid].initialize( lhs, m_inpoel, m_coord, inbox, elemblockid,
1558 [ + - ]: 66 : ue, t0, esuel.size()/4 );
1559 : : }
1560 : :
1561 : : // Transfer initial conditions from cells to nodes
1562 [ + + ]: 34172 : for (std::size_t p=0; p<npoin; ++p) { // for all mesh nodes on this chare
1563 [ + - ][ - - ]: 34106 : std::vector< tk::real > up( nprop, 0.0 );
1564 : : tk::real vol = 0.0;
1565 [ + + ]: 533606 : for (auto e : tk::Around(esup,p)) { // for all cells around node p
1566 : : // compute nodal volume: every element contributes their volume / 4
1567 : 499500 : vol += geoElem(e,0) / 4.0;
1568 : : // sum cell value to node weighed by cell volume / 4
1569 [ + + ]: 1209240 : for (std::size_t c=0; c<nprop; ++c)
1570 [ + - ]: 1419480 : up[c] += ue[e][c] * geoElem(e,0) / 4.0;
1571 : : }
1572 : : // store nodal value
1573 [ + + ]: 81412 : for (std::size_t c=0; c<nprop; ++c) u(p,c) = up[c] / vol;
1574 : : }
1575 : :
1576 [ - - ][ - - ]: 0 : } else Throw( "Scheme centring not handled for nodal initialization" );
[ - - ][ - - ]
[ - - ][ - - ]
1577 : :
1578 : : Assert( u.nunk() == m_coord[0].size(), "Size mismatch" );
1579 : : Assert( u.nprop() == nprop, "Size mismatch" );
1580 : :
1581 : 165 : return u;
1582 : : }
1583 : :
1584 : : void
1585 : 177 : Refiner::updateMesh()
1586 : : // *****************************************************************************
1587 : : // Update old mesh after refinement
1588 : : // *****************************************************************************
1589 : : {
1590 : : // Get refined mesh connectivity
1591 : 177 : const auto& refinpoel = m_refiner.tet_store.get_active_inpoel();
1592 : : Assert( refinpoel.size()%4 == 0, "Inconsistent refined mesh connectivity" );
1593 : :
1594 : : // Generate unique node lists of old and refined mesh using local ids
1595 : 177 : auto rinpoel = m_inpoel;
1596 [ + - ]: 177 : tk::remap( rinpoel, m_rid );
1597 [ + - ][ - - ]: 177 : std::unordered_set< std::size_t > old( begin(rinpoel), end(rinpoel) );
1598 : 177 : std::unordered_set< std::size_t > ref( begin(refinpoel), end(refinpoel) );
1599 : :
1600 : : // Augment refiner id -> local node id map with newly added nodes
1601 : : std::size_t l = m_lref.size();
1602 [ + + ][ + + ]: 205895 : for (auto r : ref) if (old.find(r) == end(old)) m_lref[r] = l++;
[ + - ]
1603 : :
1604 : : // Get nodal communication map from Discretization worker
1605 : 177 : if ( m_mode == RefMode::DTREF ||
1606 [ + + ]: 177 : m_mode == RefMode::OUTREF ||
1607 : : m_mode == RefMode::OUTDEREF ) {
1608 : : m_nodeCommMap =
1609 [ + - ]: 132 : m_scheme[m_meshid].disc()[thisIndex].ckLocal()->NodeCommMap();
1610 : : }
1611 : :
1612 : : // Update mesh and solution after refinement
1613 [ + - ]: 177 : newVolMesh( old, ref );
1614 : :
1615 : : // Update mesh connectivity from refiner lib, remapping refiner to local ids
1616 [ + - ][ + - ]: 177 : m_inpoel = m_refiner.tet_store.get_active_inpoel();
1617 [ + - ]: 177 : tk::remap( m_inpoel, m_lref );
1618 : :
1619 : : // Update mesh connectivity with new global node ids
1620 [ + - ]: 177 : m_ginpoel = m_inpoel;
1621 : : Assert( tk::uniquecopy(m_ginpoel).size() == m_coord[0].size(),
1622 : : "Size mismatch" );
1623 [ + - ]: 177 : tk::remap( m_ginpoel, m_gid );
1624 : :
1625 : : // Update boundary face and node information
1626 [ + - ]: 177 : newBndMesh( ref );
1627 : :
1628 : : // Augment node communication map with newly added nodes on chare-boundary
1629 [ + + ]: 177 : if (m_mode == RefMode::DTREF || m_mode == RefMode::OUTREF) {
1630 [ + + ]: 105 : for (const auto& [ neighborchare, edges ] : m_remoteEdges) {
1631 : : auto& nodes = tk::ref_find( m_nodeCommMap, neighborchare );
1632 [ + + ]: 24192 : for (const auto& e : edges) {
1633 : : // If parent nodes were part of the node communication map for chare
1634 [ + + ][ + + ]: 28404 : if (nodes.find(e[0]) != end(nodes) && nodes.find(e[1]) != end(nodes)) {
1635 : : // Add new node if local id was generated for it
1636 : 2418 : auto n = Hash<2>()( e );
1637 [ + + ]: 4836 : if (m_lid.find(n) != end(m_lid)) nodes.insert( n );
1638 : : }
1639 : : }
1640 : : }
1641 : : }
1642 : :
1643 : : // Ensure valid mesh after refinement
1644 : : Assert( tk::positiveJacobians( m_inpoel, m_coord ),
1645 : : "Refined mesh cell Jacobian non-positive" );
1646 : :
1647 : : Assert( tk::conforming( m_inpoel, m_coord, true, m_rid ),
1648 : : "Chare-"+std::to_string(thisIndex)+
1649 : : " mesh not conforming after updating mesh after mesh refinement" );
1650 : :
1651 : : // Perform leak test on new mesh
1652 : : Assert( !tk::leakyPartition(
1653 : : tk::genEsuelTet( m_inpoel, tk::genEsup(m_inpoel,4) ),
1654 : : m_inpoel, m_coord ),
1655 : : "Refined mesh partition leaky" );
1656 : 177 : }
1657 : :
1658 : : void
1659 : 177 : Refiner::newVolMesh( const std::unordered_set< std::size_t >& old,
1660 : : const std::unordered_set< std::size_t >& ref )
1661 : : // *****************************************************************************
1662 : : // Compute new volume mesh after mesh refinement
1663 : : //! \param[in] old Unique nodes of the old (unrefined) mesh using
1664 : : //! refiner-lib ids
1665 : : //! \param[in] ref Unique nodes of the refined mesh using refiner-lib ids
1666 : : // *****************************************************************************
1667 : : {
1668 : : const auto& x = m_coord[0];
1669 : : const auto& y = m_coord[1];
1670 : : const auto& z = m_coord[2];
1671 : :
1672 : : // Generate coordinates and ids to newly added nodes after refinement
1673 : : std::unordered_map< std::size_t, std::size_t > gid_add;
1674 : 177 : tk::destroy( m_addedNodes );
1675 : 177 : tk::destroy( m_removedNodes );
1676 : 177 : tk::destroy( m_amrNodeMap );
1677 [ + + ]: 103036 : for (auto r : ref) { // for all unique nodes of the refined mesh
1678 [ + + ]: 102859 : if (old.find(r) == end(old)) { // if node is newly added
1679 : : // get (local) parent ids of newly added node
1680 : 73546 : auto p = m_refiner.node_connectivity.get( r );
1681 : : Assert(p[0] != p[1], "Node without parent edge in newVolMesh");
1682 : : Assert( old.find(p[0]) != end(old) && old.find(p[1]) != end(old),
1683 : : "Parent(s) not in old mesh" );
1684 : : // local parent ids
1685 : 73546 : decltype(p) lp{{tk::cref_find(m_lref,p[0]), tk::cref_find(m_lref,p[1])}};
1686 : : // global parent ids
1687 : 73546 : decltype(p) gp{{m_gid[lp[0]], m_gid[lp[1]]}};
1688 : : // generate new global ID for newly added node
1689 : 73546 : auto g = Hash<2>()( gp );
1690 : :
1691 : : // if node added by AMR lib has not yet been added to Refiner's new mesh
1692 [ + - ]: 73546 : if (m_coordmap.find(g) == end(m_coordmap)) {
1693 : : Assert( g >= old.size(), "Hashed id overwriting old id" );
1694 : : Assert( m_lid.find(g) == end(m_lid),
1695 : : "Overwriting entry global->local node ID map" );
1696 [ + - ]: 73546 : auto l = tk::cref_find( m_lref, r );
1697 : : // store newly added node id and their parent ids (local ids)
1698 : 73546 : m_addedNodes[r] = lp; // key = r for later update to local
1699 : : // assign new node to refiner->global map
1700 [ + - ]: 73546 : gid_add[r] = g; // key = r for later search
1701 : : // assign new node to global->local map
1702 [ + - ]: 73546 : m_lid[g] = l;
1703 : : // generate and store coordinates for newly added node
1704 [ + - ]: 73546 : m_coordmap.insert( {g, {{ (x[lp[0]] + x[lp[1]])/2.0,
1705 : 73546 : (y[lp[0]] + y[lp[1]])/2.0,
1706 [ + - ]: 73546 : (z[lp[0]] + z[lp[1]])/2.0 }} } );
1707 : : }
1708 : : }
1709 : : }
1710 : 177 : tk::destroy( m_coord );
1711 : :
1712 : : // generate a node map based on oldnodes+addednodes
1713 [ + - ]: 177 : std::vector< size_t > nodeVec(m_coordmap.size());
1714 [ + + ]: 152950 : for (size_t j=0; j<nodeVec.size(); ++j) {
1715 : 152773 : nodeVec[j] = j;
1716 : : }
1717 : :
1718 : : // Remove coordinates and ids of removed nodes due to derefinement
1719 : : std::unordered_map< std::size_t, std::size_t > gid_rem;
1720 [ + + ]: 79404 : for (auto o : old) { // for all unique nodes of the old mesh
1721 [ + + ]: 79227 : if (ref.find(o) == end(ref)) { // if node is no longer in new mesh
1722 : 49914 : auto l = tk::cref_find( m_lref, o );
1723 [ + - ]: 49914 : auto g = m_gid[l];
1724 : : // store local-ids of removed nodes
1725 : : m_removedNodes.insert(l);
1726 : : // remove derefined nodes from node comm map
1727 [ + + ]: 74034 : for (auto& [neighborchare, sharednodes] : m_nodeCommMap) {
1728 [ - + ]: 24120 : if (sharednodes.find(g) != sharednodes.end()) {
1729 : : sharednodes.erase(g);
1730 : : }
1731 : : }
1732 [ + - ]: 49914 : gid_rem[l] = g;
1733 : 49914 : m_lid.erase( g );
1734 : : m_coordmap.erase( g );
1735 : : }
1736 : : }
1737 : :
1738 : : // update the node map by removing the derefined nodes
1739 [ - + ][ - - ]: 177 : if (m_mode == RefMode::DTREF && m_removedNodes.size() > 0) {
1740 : : // remove derefined nodes
1741 : : size_t remCount = 0;
1742 : 0 : size_t origSize = nodeVec.size();
1743 [ - - ]: 0 : for (size_t j=0; j<origSize; ++j) {
1744 : 0 : auto nd = nodeVec[j-remCount];
1745 : :
1746 : : bool no_change = false;
1747 : : size_t nodeidx = 0;
1748 [ - - ]: 0 : for (const auto& rn : m_removedNodes) {
1749 [ - - ]: 0 : if (nd < *m_removedNodes.cbegin()) {
1750 : : no_change = true;
1751 : : break;
1752 : : }
1753 [ - - ]: 0 : else if (nd <= rn) {
1754 : : nodeidx = rn;
1755 : : break;
1756 : : }
1757 : : }
1758 : :
1759 : : // if node is out-or-range of removed nodes list, continue with next entry
1760 [ - - ]: 0 : if (no_change)
1761 : 0 : continue;
1762 : : // if not is within range of removed nodes list, erase node appropriately
1763 [ - - ]: 0 : else if (nodeidx == nd) {
1764 : : //! Difference type for iterator/pointer arithmetics
1765 : : using diff_type = std::vector< std::size_t >::difference_type;
1766 : : nodeVec.erase(nodeVec.begin()+static_cast< diff_type >(j-remCount));
1767 : 0 : ++remCount;
1768 : : }
1769 : : }
1770 : :
1771 : : Assert(remCount == m_removedNodes.size(), "Incorrect number of nodes removed "
1772 : : "from node map.");
1773 : : }
1774 : :
1775 : : // invert node vector to get node map
1776 [ + + ]: 152950 : for (size_t i=0; i<nodeVec.size(); ++i) {
1777 : 152773 : m_amrNodeMap[nodeVec[i]] = i;
1778 : : }
1779 : :
1780 : : //// Save previous states of refiner-local node id maps before update
1781 : : //m_oldrid = m_rid;
1782 : : //m_oldlref = m_lref;
1783 : :
1784 : : // Generate new node id maps for nodes kept
1785 : 177 : tk::destroy( m_lref );
1786 [ + - ][ + - ]: 177 : std::vector< std::size_t > rid( ref.size() );
1787 [ + - ][ - - ]: 177 : std::vector< std::size_t > gid( ref.size() );
1788 : 177 : std::size_t l = 0; // will generate new local node id
1789 [ + + ]: 79404 : for (std::size_t i=0; i<m_gid.size(); ++i) {
1790 [ + + ]: 79227 : if (gid_rem.find(i) == end(gid_rem)) {
1791 : 29313 : gid[l] = m_gid[i];
1792 : 29313 : rid[l] = m_rid[i];
1793 [ + - ]: 29313 : m_lref[ m_rid[i] ] = l;
1794 : 29313 : ++l;
1795 : : }
1796 : : }
1797 : : // Add newly added nodes due to refinement to node id maps
1798 [ - - ]: 177 : decltype(m_addedNodes) addedNodes( m_addedNodes.size() );
1799 [ + + ]: 73723 : for (const auto& n : gid_add) {
1800 : 73546 : auto r = n.first;
1801 : 73546 : auto g = n.second;
1802 [ + - ]: 73546 : gid[l] = g;
1803 : 73546 : rid[l] = r;
1804 : : Assert(m_lref.find(r) == m_lref.end(), "Overwriting lref");
1805 [ + - ]: 73546 : m_lref[r] = l;
1806 : : auto it = m_addedNodes.find( r );
1807 : : Assert( it != end(m_addedNodes), "Cannot find added node" );
1808 [ + - ]: 73546 : addedNodes[l] = std::move(it->second);
1809 [ + - ][ + - ]: 73546 : addedNodes.at(l)[0] = m_amrNodeMap[addedNodes.at(l)[0]];
1810 [ + - ]: 73546 : addedNodes.at(l)[1] = m_amrNodeMap[addedNodes.at(l)[1]];
1811 : 73546 : ++l;
1812 : : }
1813 : : Assert( m_lref.size() == ref.size(), "Size mismatch" );
1814 : : //for (auto r : ref) {
1815 : : // Assert(m_lref.find(r) != m_lref.end(), "Node missing in lref");
1816 : : //}
1817 : : //const auto& int_list = m_refiner.tet_store.intermediate_list;
1818 : : //for (auto in : int_list) {
1819 : : // Assert(m_lref.find(in) != m_lref.end(), "Interm node missing in lref: "
1820 : : // + std::to_string(in));
1821 : : //}
1822 : : m_rid = std::move( rid );
1823 : : Assert( m_rid.size() == ref.size(), "Size mismatch" );
1824 : : m_addedNodes = std::move( addedNodes );
1825 : :
1826 : : // Update node coordinates, ids, and id maps
1827 : : auto& rx = m_coord[0];
1828 : : auto& ry = m_coord[1];
1829 : : auto& rz = m_coord[2];
1830 [ + - ]: 177 : rx.resize( ref.size() );
1831 [ + - ]: 177 : ry.resize( ref.size() );
1832 [ + - ]: 177 : rz.resize( ref.size() );
1833 [ + + ]: 103036 : for (std::size_t i=0; i<gid.size(); ++i) {
1834 : 205718 : tk::ref_find( m_lid, gid[i] ) = i;
1835 : : const auto& c = tk::cref_find( m_coordmap, gid[i] );
1836 : 102859 : rx[i] = c[0];
1837 : 102859 : ry[i] = c[1];
1838 : 102859 : rz[i] = c[2];
1839 : : }
1840 [ + - ]: 177 : m_gid = std::move( gid );
1841 : : Assert( m_gid.size() == m_lid.size() && m_gid.size() == ref.size(),
1842 : : "Size mismatch" );
1843 : 177 : }
1844 : :
1845 : : std::unordered_set< std::size_t >
1846 : 3743942 : Refiner::ancestors( std::size_t n )
1847 : : // *****************************************************************************
1848 : : // Find the oldest parents of a mesh node in the AMR hierarchy
1849 : : //! \param[in] n Local node id whose ancestors to search
1850 : : //! \return Parents of local node id from the coarsest (original) mesh
1851 : : // *****************************************************************************
1852 : : {
1853 : 3743942 : auto d = m_refiner.node_connectivity.get( m_rid[n] );
1854 : 3743942 : decltype(d) p{{ tk::cref_find( m_lref, d[0] ),
1855 [ + + ]: 7487884 : tk::cref_find( m_lref, d[1] ) }};
1856 : :
1857 : : std::unordered_set< std::size_t > s;
1858 : :
1859 [ + + ]: 3743942 : if (p != AMR::node_pair_t{{n,n}}) {
1860 [ + - ]: 1072675 : auto q = ancestors( p[0] );
1861 : : s.insert( begin(q), end(q) );
1862 [ + - ]: 1072675 : auto r = ancestors( p[1] );
1863 : : s.insert( begin(r), end(r) );
1864 : : } else {
1865 : : s.insert( begin(p), end(p) );
1866 : : }
1867 : :
1868 : 3743942 : return s;
1869 : : }
1870 : :
1871 : : Refiner::BndFaceData
1872 : 177 : Refiner::boundary()
1873 : : // *****************************************************************************
1874 : : // Generate boundary data structures used to update refined/derefined boundary
1875 : : // faces and nodes of side sets
1876 : : //! \return A tuple of boundary face data
1877 : : //! \details The output of this function is used to regenerate physical boundary
1878 : : //! face and node data structures after refinement, see updateBndData().
1879 : : // *****************************************************************************
1880 : : {
1881 : : // Generate the inverse of AMR's tet store.
1882 : : std::unordered_map< Tet, std::size_t, Hash<4>, Eq<4> > invtets;
1883 [ + + ]: 439679 : for (const auto& [key, tet] : m_refiner.tet_store.tets)
1884 [ + - ]: 439502 : invtets[ tet ] = key;
1885 : :
1886 : : //std::cout << thisIndex << " invt: " << invtets.size() << '\n';
1887 : : //std::cout << thisIndex << " active inpoel size: " << m_refiner.tet_store.get_active_inpoel().size() << '\n';
1888 : : //std::cout << thisIndex << " marked derefinement size: " << m_refiner.tet_store.marked_derefinements.size() << '\n';
1889 : :
1890 : : // Generate data structure pcFaceTets for the new (post-AMR) mesh:
1891 : : // pcFaceTets is a map that associates all triangle boundary faces (physical
1892 : : // and chare) to the id of the tet adjacent to the said face.
1893 : : // Key: Face-nodes' global id; Value: tet-id.
1894 : : std::unordered_map< Face, std::size_t, Hash<3>, Eq<3> > pcFaceTets;
1895 [ + - ][ + - ]: 354 : auto esuel = tk::genEsuelTet( m_inpoel, tk::genEsup(m_inpoel,4) );
1896 [ + + ]: 392662 : for (std::size_t e=0; e<esuel.size()/4; ++e) {
1897 : 392485 : auto m = e*4;
1898 [ + + ]: 1962425 : for (std::size_t f=0; f<4; ++f) {
1899 [ + + ]: 1569940 : if (esuel[m+f] == -1) { // if a face does not have an adjacent tet
1900 : 139186 : Face b{{ m_ginpoel[ m+tk::lpofa[f][0] ],
1901 : 139186 : m_ginpoel[ m+tk::lpofa[f][1] ],
1902 : 139186 : m_ginpoel[ m+tk::lpofa[f][2] ] }};
1903 : : Assert( m_inpoel[m+0] < m_rid.size() &&
1904 : : m_inpoel[m+1] < m_rid.size() &&
1905 : : m_inpoel[m+2] < m_rid.size() &&
1906 : : m_inpoel[m+3] < m_rid.size(), "Indexing out of rid" );
1907 : 139186 : Tet t{{ m_rid[ m_inpoel[m+0] ], m_rid[ m_inpoel[m+1] ],
1908 : 139186 : m_rid[ m_inpoel[m+2] ], m_rid[ m_inpoel[m+3] ] }};
1909 : : //Tet t{{ m_inpoel[m+0], m_inpoel[m+1],
1910 : : // m_inpoel[m+2], m_inpoel[m+3] }};
1911 : : // associate tet id to adjacent (physical or chare) boundary face
1912 : : auto i = invtets.find( t );
1913 : : Assert(m_refiner.tet_store.is_active(i->second),
1914 : : "Inactive element while regenerating boundary data");
1915 [ + - ]: 139186 : if (i != end(invtets)) {
1916 : : //std::cout << "refacetets: " <<
1917 : : // b[0] << "-" << b[1] << "-" << b[2] << std::endl;
1918 [ + - ]: 139186 : pcFaceTets[ b ] = i->second;
1919 : : } else {
1920 [ - - ][ - - ]: 0 : Throw("Active element not found in tet_store");
[ - - ][ - - ]
[ - - ][ - - ]
1921 : : }
1922 : : }
1923 : : }
1924 : : }
1925 : :
1926 : : // Generate child->parent tet and id maps after refinement/derefinement step
1927 : : // tk::destroy( m_oldparent );
1928 : : m_addedTets.clear();
1929 : : std::size_t p = 0;
1930 : : std::size_t c = 0;
1931 : : const auto& tet_store = m_refiner.tet_store;
1932 [ + + ]: 439679 : for (const auto& t : tet_store.tets) {
1933 : : // query number of children of tet
1934 [ + - ]: 439502 : auto nc = tet_store.data( t.first ).children.size();
1935 [ + + ]: 813096 : for (decltype(nc) i=0; i<nc; ++i ) { // for all child tets
1936 : : // get child tet id
1937 [ + - ]: 373594 : auto childtet = tet_store.get_child_id( t.first, i );
1938 : : auto ct = tet_store.tets.find( childtet );
1939 : : Assert(ct != tet_store.tets.end(), "Child not found in tet store");
1940 : : // //auto cA = tk::cref_find( m_lref, ct->second[0] );
1941 : : // //auto cB = tk::cref_find( m_lref, ct->second[1] );
1942 : : // //auto cC = tk::cref_find( m_lref, ct->second[2] );
1943 : : // //auto cD = tk::cref_find( m_lref, ct->second[3] );
1944 : : // // get nodes of parent tet
1945 : : // //auto pA = tk::cref_find( m_lref, t.second[0] );
1946 : : // //auto pB = tk::cref_find( m_lref, t.second[1] );
1947 : : // //auto pC = tk::cref_find( m_lref, t.second[2] );
1948 : : // //auto pD = tk::cref_find( m_lref, t.second[3] );
1949 : : // // assign parent tet to child tet
1950 : : // //m_oldparent[ {{cA,cB,cC,cD}} ] = {{pA,pB,pC,pD}};
1951 : : // m_oldparent[ ct->second ] = t.second; //{{pA,pB,pC,pD}};
1952 [ + + ]: 373594 : if (m_oldTets.find(ct->second) == end(m_oldTets)) {
1953 : : // TODO: the following code can assign negative ids to newly added tets.
1954 : : // This needs to be corrected before applying to cell-based schemes.
1955 : : //Assert((p-m_oldntets) > 0, "Negative id assigned to added tet");
1956 [ + - ][ - - ]: 350030 : m_addedTets[ c++ ] = p - m_oldntets;
1957 : : }
1958 : : }
1959 : 439502 : ++p;
1960 : : }
1961 : :
1962 : : //std::cout << thisIndex << " added: " << m_addedTets.size() << '\n';
1963 : : //std::cout << thisIndex << " parent: " << m_oldparent.size() << '\n';
1964 : : //std::cout << thisIndex << " pcret: " << pcFaceTets.size() << '\n';
1965 : :
1966 : : //for (std::size_t f=0; f<m_triinpoel.size()/3; ++f) {
1967 : : // std::cout << "triinpoel: " <<
1968 : : // m_triinpoel[f*3+0] << "-" << m_triinpoel[f*3+1] << "-" <<
1969 : : // m_triinpoel[f*3+2] << std::endl;
1970 : : //}
1971 : :
1972 : 177 : return pcFaceTets;
1973 : : }
1974 : :
1975 : : void
1976 : 177 : Refiner::newBndMesh( const std::unordered_set< std::size_t >& ref )
1977 : : // *****************************************************************************
1978 : : // Update boundary data structures after mesh refinement
1979 : : //! \param[in] ref Unique nodes of the refined mesh using refiner-lib ids
1980 : : // *****************************************************************************
1981 : : {
1982 : : // Generate boundary face data structures used to regenerate boundary face
1983 : : // and node data after mesh refinement
1984 : 177 : auto pcFaceTets = boundary();
1985 : :
1986 : : // Regerate boundary faces and nodes after AMR step
1987 [ + - ]: 177 : updateBndData( ref, pcFaceTets );
1988 : 177 : }
1989 : :
1990 : : void
1991 : 177 : Refiner::updateBndData(
1992 : : [[maybe_unused]] const std::unordered_set< std::size_t >& ref,
1993 : : const BndFaceData& pcFaceTets )
1994 : : // *****************************************************************************
1995 : : // Regenerate boundary faces and nodes after AMR step
1996 : : //! \param[in] ref Unique nodes of the refined mesh using refiner-lib ids
1997 : : //! \param[in] pcFaceTets Boundary face data
1998 : : // *****************************************************************************
1999 : : {
2000 : : // storage for boundary faces associated to side-set IDs of the refined mesh
2001 : 177 : tk::destroy( m_bface );
2002 : : // storage for boundary faces-node connectivity of the refined mesh
2003 : : tk::destroy( m_triinpoel );
2004 : : // storage for boundary nodes associated to side-set IDs of the refined mesh
2005 : 177 : tk::destroy( m_bnode );
2006 : :
2007 : : // face id counter
2008 : 177 : std::size_t facecnt = 0;
2009 : : // will collect unique faces added for each side set
2010 : : std::unordered_map< int, FaceSet > bf;
2011 : :
2012 : : // Lambda to associate a boundary face and connectivity to a side set.
2013 : : // Argument 's' is the list of faces (ids) to add the new face to. Argument
2014 : : // 'ss' is the side set id to which the face is added. Argument 'f' is the
2015 : : // triangle face connectivity to add.
2016 : 103510 : auto addBndFace = [&]( std::vector< std::size_t >& s, int ss, const Face& f )
2017 : : {
2018 : : // only add face if it has not yet been aded to this side set
2019 [ + - ]: 207020 : if (bf[ ss ].insert( f ).second) {
2020 : 103510 : s.push_back( facecnt++ );
2021 : 103510 : m_triinpoel.insert( end(m_triinpoel), begin(f), end(f) );
2022 : : Assert(m_triinpoel.size()/3 == facecnt, "Incorrect size of triinpoel");
2023 : : }
2024 : 103687 : };
2025 : :
2026 : : // Lambda to search the parents in the coarsest mesh of a mesh node and if
2027 : : // found, add its global id to boundary node lists associated to the side
2028 : : // set(s) of its parents. Argument 'n' is the local id of the mesh node id
2029 : : // whose parents to search.
2030 : 1181034 : auto addBndNodes = [&]( std::size_t n ){
2031 : 1181034 : auto a = ancestors( n ); // find parents of n in coarse mesh
2032 [ + + ]: 1181034 : if (a.size() == 1) {
2033 : : // node was part of the coarse mesh
2034 : : Assert(*a.cbegin() == n, "Single ancestor not self");
2035 [ + - ]: 434862 : auto ss = keys( m_coarseBndNodes, m_gid[*a.cbegin()] );
2036 [ + + ]: 442398 : for (auto s : ss)
2037 [ + - ][ + - ]: 7536 : m_bnode[ s ].push_back( m_gid[n] );
2038 [ + - ]: 746172 : } else if (a.size() == 2) {
2039 : : // node was added to an edge of a coarse face
2040 [ + - ]: 746172 : std::vector< std::size_t > p( begin(a), end(a) );
2041 [ + - ]: 746172 : auto ss1 = keys( m_coarseBndNodes, m_gid[p[0]] );
2042 [ + - ]: 746172 : auto ss2 = keys( m_coarseBndNodes, m_gid[p[1]] );
2043 [ + + ]: 765996 : for (auto s : ss1) {
2044 : : // only add 'n' to bnode if all parent nodes are in same side set, else
2045 : : // 'n' is not a boundary node
2046 [ + + ]: 19824 : if (ss2.find(s) != end(ss2)) {
2047 [ + - ][ + - ]: 17280 : m_bnode[ s ].push_back( m_gid[n] );
2048 : : }
2049 : : }
2050 [ - - ]: 0 : } else if (a.size() == 3) {
2051 : : // node was added inside of a coarse face
2052 [ - - ]: 0 : std::vector< std::size_t > p( begin(a), end(a) );
2053 [ - - ]: 0 : auto ss1 = keys( m_coarseBndNodes, m_gid[p[0]] );
2054 [ - - ]: 0 : auto ss2 = keys( m_coarseBndNodes, m_gid[p[1]] );
2055 [ - - ]: 0 : auto ss3 = keys( m_coarseBndNodes, m_gid[p[2]] );
2056 [ - - ]: 0 : for (auto s : ss1) {
2057 : : // only add 'n' to bnode if all parent nodes are in same side set, else
2058 : : // 'n' is not a boundary node
2059 [ - - ][ - - ]: 0 : if (ss2.find(s) != end(ss2) && ss3.find(s) != end(ss3)) {
2060 [ - - ][ - - ]: 0 : m_bnode[ s ].push_back( m_gid[n] );
2061 : : }
2062 : : }
2063 : : }
2064 : 1181211 : };
2065 : :
2066 : : // Regenerate boundary faces for new mesh along side sets. For all faces
2067 : : // associated to side sets, we find the ancestors (parents of nodes in the
2068 : : // original/coarsest mesh) of the nodes comprising the physical and chare
2069 : : // boundary faces of the new mesh.
2070 : : //bool faceNoSs = false;
2071 : : // for all P/C faces in current (post-AMR) mesh
2072 [ + + ]: 139363 : for (const auto& [ face, tetid ] : pcFaceTets) {
2073 : : // find ancestors of face
2074 : : std::unordered_set< std::size_t > ans;
2075 [ + + ]: 556744 : for (std::size_t i=0; i<3; ++i) {
2076 [ + - ]: 835116 : auto ai = ancestors(tk::cref_find(m_lid, face[i]));
2077 : : ans.insert(ai.begin(), ai.end());
2078 : : }
2079 : : Assert(ans.size() == 3, "Incorrect number of ancestors in refined face");
2080 : : Face af;
2081 : : std::size_t i = 0;
2082 [ + + ]: 556744 : for (auto ai:ans) {
2083 : 417558 : af[i] = m_gid[ai];
2084 : 417558 : ++i;
2085 : : }
2086 : : // for all boundary faces in original mesh
2087 : : //std::size_t fss = 0;
2088 [ + + ]: 532864 : for (const auto& [ss, cfaceset] : m_coarseBndFaces) {
2089 [ + + ]: 393678 : if (cfaceset.find(af) != cfaceset.end()) {
2090 [ + - ][ + - ]: 103510 : addBndFace(m_bface[ss], ss, face);
2091 : : //++fss;
2092 : : }
2093 [ + + ]: 1574712 : for (auto j : face) {
2094 [ + - ]: 2362068 : addBndNodes(tk::cref_find(m_lid, j));
2095 : : }
2096 : : }
2097 : : //if (fss==0) {
2098 : : // std::cout << "Face added to no/multiple side sets; " << fss << std::endl;
2099 : : // faceNoSs = true;
2100 : : //}
2101 : : }
2102 : :
2103 : : // Commented code below, since diagcg can work without sideset/bcs
2104 : : //Assert(!faceNoSs, "Face/s added to incorrect number of side sets");
2105 : :
2106 : : // Make boundary node IDs unique for each physical boundary (side set)
2107 [ + + ][ + - ]: 193 : for (auto& s : m_bnode) tk::unique( s.second );
2108 : :
2109 : : //for (const auto& [ setid, faceids ] : m_bface) {
2110 : : // std::cout << "sset: " << setid << std::endl;
2111 : : // for (auto f : faceids) {
2112 : : // Assert(f<m_triinpoel.size()/3, "Out of bounds access into triinpoel");
2113 : : // std::cout << "new bndfaces: " <<
2114 : : // m_triinpoel[f*3+0] << "-" << m_triinpoel[f*3+1] << "-" <<
2115 : : // m_triinpoel[f*3+2] << std::endl;
2116 : : // }
2117 : : //}
2118 : :
2119 : : //for (std::size_t f=0; f<m_triinpoel.size()/3; ++f) {
2120 : : // std::cout << "new triinpoel: " <<
2121 : : // m_triinpoel[f*3+0] << "-" << m_triinpoel[f*3+1] << "-" <<
2122 : : // m_triinpoel[f*3+2] << std::endl;
2123 : : //}
2124 : :
2125 : : //std::cout << thisIndex << " bf: " << tk::sumvalsize( m_bface ) << '\n';
2126 : :
2127 : : //std::cout << thisIndex << " bn: " << tk::sumvalsize( m_bnode ) << '\n';
2128 : :
2129 : : // Perform leak-test on boundary face data just updated (only in DEBUG)
2130 : : Assert( bndIntegral(), "Partial boundary integral" );
2131 : 177 : }
2132 : :
2133 : : bool
2134 : 0 : Refiner::bndIntegral()
2135 : : // *****************************************************************************
2136 : : // Compute partial boundary surface integral and sum across all chares
2137 : : //! \return true so we don't trigger assert in client code
2138 : : //! \details This function computes a partial surface integral over the boundary
2139 : : //! of the faces of this mesh partition then sends its contribution to perform
2140 : : //! the integral acorss the total problem boundary. After the global sum a
2141 : : //! non-zero vector result indicates a leak, e.g., a hole in the boundary
2142 : : //! which indicates an error in the boundary face data structures used to
2143 : : //! compute the partial surface integrals.
2144 : : // *****************************************************************************
2145 : : {
2146 : : const auto& x = m_coord[0];
2147 : : const auto& y = m_coord[1];
2148 : : const auto& z = m_coord[2];
2149 : :
2150 : 0 : std::vector< tk::real > s{{ 0.0, 0.0, 0.0 }};
2151 : :
2152 [ - - ]: 0 : for (const auto& [ setid, faceids ] : m_bface) {
2153 [ - - ]: 0 : for (auto f : faceids) {
2154 : 0 : auto A = tk::cref_find( m_lid, m_triinpoel[f*3+0] );
2155 : 0 : auto B = tk::cref_find( m_lid, m_triinpoel[f*3+1] );
2156 : 0 : auto C = tk::cref_find( m_lid, m_triinpoel[f*3+2] );
2157 : : // Compute geometry data for face
2158 : 0 : auto geoface = tk::geoFaceTri( {{x[A], x[B], x[C]}},
2159 [ - - ]: 0 : {{y[A], y[B], y[C]}},
2160 [ - - ]: 0 : {{z[A], z[B], z[C]}} );
2161 : : // Sum up face area * face unit-normal
2162 : 0 : s[0] += geoface(0,0) * geoface(0,1);
2163 : 0 : s[1] += geoface(0,0) * geoface(0,2);
2164 : 0 : s[2] += geoface(0,0) * geoface(0,3);
2165 : : }
2166 : : }
2167 : :
2168 [ - - ]: 0 : s.push_back( -1.0 ); // negative: no call-back after reduction
2169 [ - - ][ - - ]: 0 : s.push_back( static_cast< tk::real >( m_meshid ) );
2170 : :
2171 : : // Send contribution to host summing partial surface integrals
2172 [ - - ]: 0 : contribute( s, CkReduction::sum_double, m_cbr.get< tag::bndint >() );
2173 : :
2174 : 0 : return true; // don't trigger the assert in client code
2175 : : }
2176 : :
2177 : : #include "NoWarning/refiner.def.h"
|
