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