Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/Discretization.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \details Data and functionality common to all discretization schemes
9 : : \see Discretization.h and Discretization.C for more info.
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "Tags.hpp"
14 : : #include "Reorder.hpp"
15 : : #include "Vector.hpp"
16 : : #include "DerivedData.hpp"
17 : : #include "Discretization.hpp"
18 : : #include "MeshWriter.hpp"
19 : : #include "DiagWriter.hpp"
20 : : #include "Inciter/InputDeck/InputDeck.hpp"
21 : : #include "Inciter/Options/Scheme.hpp"
22 : : #include "Print.hpp"
23 : : #include "Around.hpp"
24 : : #include "QuinoaBuildConfig.hpp"
25 : : #include "ConjugateGradients.hpp"
26 : : #include "ALE.hpp"
27 : :
28 : : #include "M2MTransfer.hpp"
29 : :
30 : : namespace inciter {
31 : :
32 : : static CkReduction::reducerType PDFMerger;
33 : : extern ctr::InputDeck g_inputdeck;
34 : : extern ctr::InputDeck g_inputdeck_defaults;
35 : :
36 : : } // inciter::
37 : :
38 : : using inciter::Discretization;
39 : :
40 : 1943 : Discretization::Discretization(
41 : : std::size_t meshid,
42 : : const std::vector< CProxy_Discretization >& disc,
43 : : const CProxy_ALE& aleproxy,
44 : : const tk::CProxy_ConjugateGradients& conjugategradientsproxy,
45 : : const CProxy_Transporter& transporter,
46 : : const tk::CProxy_MeshWriter& meshwriter,
47 : : const tk::UnsMesh::CoordMap& coordmap,
48 : : const tk::UnsMesh::Chunk& el,
49 : : const tk::CommMaps& msum,
50 : : const std::map< int, std::vector< std::size_t > >& bface,
51 : : const std::vector< std::size_t >& triinpoel,
52 : : const std::unordered_map< std::size_t, std::set< std::size_t > >& elemblockid,
53 : 1943 : int nc ) :
54 : : m_meshid( meshid ),
55 : 1943 : m_transfer( g_inputdeck.get< tag::transfer >() ),
56 : : m_disc( disc ),
57 : : m_nchare( nc ),
58 : : m_it( 0 ),
59 : : m_itr( 0 ),
60 : : m_itf( 0 ),
61 : : m_initial( 1 ),
62 : 3886 : m_t( g_inputdeck.get< tag::t0 >() ),
63 : 3886 : m_lastDumpTime( -std::numeric_limits< tk::real >::max() ),
64 : : m_physFieldFloor( 0.0 ),
65 : : m_physHistFloor( 0.0 ),
66 : : m_rangeFieldFloor( 0.0 ),
67 : : m_rangeHistFloor( 0.0 ),
68 : 3886 : m_dt( g_inputdeck.get< tag::dt >() ),
69 : 1943 : m_dtn( m_dt ),
70 : : m_nvol( 0 ),
71 : : m_nxfer( 0 ),
72 : : m_ale( aleproxy ),
73 : : m_transporter( transporter ),
74 : : m_meshwriter( meshwriter ),
75 : : m_el( el ), // fills m_inpoel, m_gid, m_lid
76 : : m_coord( setCoord( coordmap ) ),
77 : 1943 : m_coordn( m_coord ),
78 : : m_nodeCommMap(),
79 : : m_edgeCommMap(),
80 : : m_meshvol( 0.0 ),
81 : 0 : m_v( m_gid.size(), 0.0 ),
82 : 0 : m_vol( m_gid.size(), 0.0 ),
83 : : m_volc(),
84 : 1943 : m_voln( m_vol ),
85 : 1943 : m_vol0( m_inpoel.size()/4, 0.0 ),
86 : : m_bid(),
87 : : m_timer(),
88 : : m_refined( 0 ),
89 : 1943 : m_prevstatus( std::chrono::high_resolution_clock::now() ),
90 : : m_nrestart( 0 ),
91 : : m_histdata(),
92 : : m_nsrc( 0 ),
93 : : m_ndst( 0 ),
94 : : m_meshvel( 0, 3 ),
95 : : m_meshvel_converged( true ),
96 : : m_bface( bface ),
97 : : m_triinpoel( triinpoel ),
98 [ + - ][ + - ]: 3886 : m_elemblockid( elemblockid )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
99 : : // *****************************************************************************
100 : : // Constructor
101 : : //! \param[in] meshid Mesh ID
102 : : //! \param[in] disc All Discretization proxies (one per mesh)
103 : : //! \param[in] aleproxy Distributed ALE proxy
104 : : //! \param[in] conjugategradientsproxy Distributed Conjugrate Gradients linear
105 : : //! solver proxy
106 : : //! \param[in] transporter Host (Transporter) proxy
107 : : //! \param[in] meshwriter Mesh writer proxy
108 : : //! \param[in] coordmap Coordinates of mesh nodes and their global IDs
109 : : //! \param[in] el Elements of the mesh chunk we operate on
110 : : //! \param[in] msum Communication maps associated to chare IDs bordering the
111 : : //! mesh chunk we operate on
112 : : //! \param[in] bface Face lists mapped to side set ids
113 : : //! \param[in] triinpoel Interconnectivity of points and boundary-faces
114 : : //! \param[in] elemblockid Local tet ids associated with mesh block ids
115 : : //! \param[in] nc Total number of Discretization chares
116 : : // *****************************************************************************
117 : : {
118 [ - + ][ - - ]: 1943 : Assert( !m_inpoel.empty(), "No elements assigned to Discretization chare" );
[ - - ][ - - ]
119 [ + - ][ - + ]: 1943 : Assert( tk::positiveJacobians( m_inpoel, m_coord ),
[ - - ][ - - ]
[ - - ]
120 : : "Jacobian in input mesh to Discretization non-positive" );
121 : : #if not defined(__INTEL_COMPILER) || defined(NDEBUG)
122 : : // The above ifdef skips running the conformity test with the intel compiler
123 : : // in debug mode only. This is necessary because in tk::conforming(), filling
124 : : // up the map can fail with some meshes (only in parallel), e.g., tube.exo,
125 : : // used by some regression tests, due to the intel compiler generating some
126 : : // garbage incorrect code - only in debug, only in parallel, only with that
127 : : // mesh.
128 [ + - ][ - + ]: 1943 : Assert( tk::conforming( m_inpoel, m_coord ),
[ - - ][ - - ]
[ - - ]
129 : : "Input mesh to Discretization not conforming" );
130 : : #endif
131 : :
132 : : // Store communication maps
133 [ + + ]: 21595 : for (const auto& [ c, maps ] : msum) {
134 [ + - ][ + - ]: 19652 : m_nodeCommMap[c] = maps.get< tag::node >();
135 [ + - ][ + - ]: 19652 : m_edgeCommMap[c] = maps.get< tag::edge >();
136 : : }
137 : :
138 : : // Get ready for computing/communicating nodal volumes
139 [ + - ]: 1943 : startvol();
140 : :
141 : : // Get chare-boundary node-id map
142 [ + - ]: 1943 : m_bid = genBid();
143 : :
144 : : // Find host elements of user-specified points where time histories are
145 : : // saved, and save the shape functions evaluated at the point locations
146 : 1943 : const auto& pt = g_inputdeck.get< tag::history_output, tag::point >();
147 [ + + ]: 1995 : for (std::size_t p=0; p<pt.size(); ++p) {
148 : : std::array< tk::real, 4 > N;
149 : 52 : const auto& l = pt[p].get< tag::coord >();
150 : 52 : const auto& id = pt[p].get< tag::id >();
151 [ + + ]: 25008 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
152 [ + - ][ + + ]: 24966 : if (tk::intet( m_coord, m_inpoel, l, e, N )) {
153 [ + - ][ + - ]: 10 : m_histdata.push_back( HistData{{ id, e, {l[0],l[1],l[2]}, N }} );
154 : 10 : break;
155 : : }
156 : : }
157 : : }
158 : :
159 : : // Insert ConjugrateGradients solver chare array element if needed
160 [ + + ]: 1943 : if (g_inputdeck.get< tag::ale, tag::ale >()) {
161 [ + - ]: 62 : m_ale[ thisIndex ].insert( conjugategradientsproxy,
162 : 31 : m_coord, m_inpoel,
163 [ + - ]: 31 : m_gid, m_lid, m_nodeCommMap );
164 : : } else {
165 [ + - ]: 1912 : m_meshvel.resize( m_gid.size() );
166 : : }
167 : :
168 : : // Register mesh with mesh-transfer lib
169 [ + + ][ - + ]: 1943 : if (m_disc.size() == 1 || m_transfer.empty()) {
[ + + ]
170 : : // skip transfer if single mesh or if not involved in coupling
171 [ + - ]: 1929 : transferInit();
172 : : } else {
173 [ + + ]: 14 : if (thisIndex == 0) {
174 [ + - ][ + - ]: 4 : exam2m::addMesh( thisProxy, m_nchare,
175 [ + - ][ + - ]: 8 : CkCallback( CkIndex_Discretization::transferInit(), thisProxy ) );
176 : : //std::cout << "Disc: " << m_meshid << " m2m::addMesh()\n";
177 : : }
178 : : }
179 : 1943 : }
180 : :
181 : : std::unordered_map< std::size_t, std::size_t >
182 : 1943 : Discretization::genBid()
183 : : // *****************************************************************************
184 : : // Generate the Bid data-structure based on the node communication-map
185 : : // *****************************************************************************
186 : : {
187 : : // Count the number of mesh nodes at which we receive data from other chares
188 : : // and compute map associating boundary-chare node ID to global node ID
189 [ + - ]: 3886 : std::vector< std::size_t > c( tk::sumvalsize( m_nodeCommMap ) );
190 : 1943 : std::size_t j = 0;
191 [ + + ][ + + ]: 148219 : for (const auto& [ch,n] : m_nodeCommMap) for (auto i : n) c[j++] = i;
192 [ + - ]: 1943 : tk::unique( c );
193 [ + - ]: 3886 : return tk::assignLid( c );
194 : : }
195 : :
196 : : void
197 : 1943 : Discretization::transferInit()
198 : : // *****************************************************************************
199 : : // Our mesh has been registered with the mesh-to-mesh transfer library (if
200 : : // coupled to other solver)
201 : : // *****************************************************************************
202 : : {
203 : : // Compute number of mesh points owned
204 : 1943 : std::size_t npoin = m_gid.size();
205 [ + + ][ + - ]: 187866 : for (auto g : m_gid) if (tk::slave(m_nodeCommMap,g,thisIndex)) --npoin;
[ + + ]
206 : :
207 : : // Tell the RTS that the Discretization chares have been created and compute
208 : : // the total number of mesh points across the distributed mesh
209 [ + - ]: 1943 : std::vector< std::size_t > meshdata{ m_meshid, npoin };
210 [ + - ]: 1943 : contribute( meshdata, CkReduction::sum_ulong,
211 [ + - ][ + - ]: 3886 : CkCallback( CkReductionTarget(Transporter,disccreated), m_transporter ) );
212 : 1943 : }
213 : :
214 : : void
215 : 43686 : Discretization::meshvelStart(
216 : : const tk::UnsMesh::Coords vel,
217 : : const std::vector< tk::real >& soundspeed,
218 : : const std::unordered_map< int,
219 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > >& bnorm,
220 : : tk::real adt,
221 : : CkCallback done ) const
222 : : // *****************************************************************************
223 : : // Start computing new mesh velocity for ALE mesh motion
224 : : //! \param[in] vel Fluid velocity at mesh nodes
225 : : //! \param[in] soundspeed Speed of sound at mesh nodes
226 : : //! \param[in] bnorm Face normals in boundary points associated to side sets
227 : : //! \param[in] adt alpha*dt of the RK time step
228 : : //! \param[in] done Function to continue with when mesh velocity has been
229 : : //! computed
230 : : // *****************************************************************************
231 : : {
232 [ + + ]: 43686 : if (g_inputdeck.get< tag::ale, tag::ale >())
233 [ + - ][ + - ]: 2222 : m_ale[ thisIndex ].ckLocal()->start( vel, soundspeed, done,
234 [ + - ][ + - ]: 1111 : m_coord, m_coordn, m_vol0, m_vol, bnorm, m_initial, m_it, m_t, adt );
235 : : else
236 : 42575 : done.send();
237 : 43686 : }
238 : :
239 : : const tk::Fields&
240 : 44304 : Discretization::meshvel() const
241 : : // *****************************************************************************
242 : : //! Query the mesh velocity
243 : : //! \return Mesh velocity
244 : : // *****************************************************************************
245 : : {
246 [ + + ]: 44304 : if (g_inputdeck.get< tag::ale, tag::ale >())
247 [ + - ]: 2235 : return m_ale[ thisIndex ].ckLocal()->meshvel();
248 : : else
249 : 42069 : return m_meshvel;
250 : : }
251 : :
252 : : void
253 : 43686 : Discretization::meshvelBnd(
254 : : const std::map< int, std::vector< std::size_t > >& bface,
255 : : const std::map< int, std::vector< std::size_t > >& bnode,
256 : : const std::vector< std::size_t >& triinpoel) const
257 : : // *****************************************************************************
258 : : // Query ALE mesh velocity boundary condition node lists and node lists at
259 : : // which ALE moves boundaries
260 : : // *****************************************************************************
261 : : {
262 [ + + ]: 43686 : if (g_inputdeck.get< tag::ale, tag::ale >())
263 [ + - ][ + - ]: 1111 : m_ale[ thisIndex ].ckLocal()->meshvelBnd( bface, bnode, triinpoel );
264 : 43686 : }
265 : :
266 : : void
267 : 43686 : Discretization::meshvelConv()
268 : : // *****************************************************************************
269 : : //! Assess and record mesh velocity linear solver convergence
270 : : // *****************************************************************************
271 : : {
272 : 43686 : auto smoother = g_inputdeck.get< tag::ale, tag::smoother >();
273 : :
274 [ + + ][ + + ]: 44120 : if (g_inputdeck.get< tag::ale, tag::ale >() &&
[ + + ]
275 [ + + ]: 434 : (smoother == ctr::MeshVelocitySmootherType::LAPLACE or
276 : : smoother == ctr::MeshVelocitySmootherType::HELMHOLTZ))
277 : : {
278 [ + - ][ + - ]: 801 : m_meshvel_converged &= m_ale[ thisIndex ].ckLocal()->converged();
279 : : }
280 : 43686 : }
281 : :
282 : : void
283 : 844 : Discretization::comfinal()
284 : : // *****************************************************************************
285 : : // Finish setting up communication maps and solution transfer callbacks
286 : : // *****************************************************************************
287 : : {
288 : : // Generate own subset of solver/mesh transfer list
289 [ + + ]: 858 : for (const auto& t : m_transfer) {
290 [ + + ][ + - ]: 14 : if (t.src == m_meshid || t.dst == m_meshid) {
291 [ + - ]: 14 : m_mytransfer.push_back( t );
292 : : }
293 : : }
294 : :
295 : : // Signal the runtime system that the workers have been created
296 [ + - ]: 844 : std::vector< std::size_t > meshdata{ /* initial = */ 1, m_meshid };
297 [ + - ]: 844 : contribute( meshdata, CkReduction::sum_ulong,
298 [ + - ][ + - ]: 1688 : CkCallback(CkReductionTarget(Transporter,comfinal), m_transporter) );
299 : 844 : }
300 : :
301 : : void
302 : 12264 : Discretization::transfer(
303 : : tk::Fields& u,
304 : : std::size_t dirn,
305 : : CkCallback cb )
306 : : // *****************************************************************************
307 : : // Start solution transfer (if coupled)
308 : : //! \param[in,out] u Solution to transfer from/to
309 : : //! \param[in] dirn Direction of solution transfer. 0: from background to
310 : : //! overset, 1: from overset to background
311 : : //! \param[in] cb Callback to call when back and forth transfers complete.
312 : : //! \details This function initiates the solution transfer (direction dependent
313 : : //! on 'dirn') between meshes. It invokes a reduction to Transporter when the
314 : : //! transfer in one direction is complete (dirn == 0), or calls back the
315 : : //! 'cb' function in Scheme when transfers both directions are complete.
316 : : //! The function relies on 'dirn' to make this decision.
317 : : // *****************************************************************************
318 : : {
319 [ + + ]: 12264 : if (m_mytransfer.empty()) { // skip transfer if not involved in coupling
320 : :
321 : 9800 : cb.send();
322 : :
323 : : } else {
324 : :
325 : 2464 : m_transfer_complete = cb;
326 : :
327 : : // determine source and destination mesh depending on direction of transfer
328 : 2464 : std::size_t fromMesh(0), toMesh(0);
329 : 4928 : CkCallback cb_xfer;
330 [ + + ]: 2464 : if (dirn == 0) {
331 : 1232 : fromMesh = m_mytransfer[m_nsrc].src;
332 : 1232 : toMesh = m_mytransfer[m_ndst].dst;
333 [ + - ][ + - ]: 1232 : cb_xfer = CkCallback( CkIndex_Discretization::to_complete(), thisProxy[thisIndex] );
[ + - ]
334 : : }
335 : : else {
336 : 1232 : fromMesh = m_mytransfer[m_nsrc].dst;
337 : 1232 : toMesh = m_mytransfer[m_ndst].src;
338 [ + - ][ + - ]: 1232 : cb_xfer = CkCallback( CkIndex_Discretization::from_complete(), thisProxy[thisIndex] );
[ + - ]
339 : : }
340 : :
341 : : // Pass source and destination meshes to mesh transfer lib (if coupled)
342 [ - + ][ - - ]: 2464 : Assert( m_nsrc < m_mytransfer.size(), "Indexing out of mytransfer[src]" );
[ - - ][ - - ]
343 [ + + ]: 2464 : if (fromMesh == m_meshid) {
344 [ + - ][ + - ]: 1232 : exam2m::setSourceTets( thisProxy, thisIndex, &m_inpoel, &m_coord, u );
345 : 1232 : ++m_nsrc;
346 : : } else {
347 : 1232 : m_nsrc = 0;
348 : : }
349 [ - + ][ - - ]: 2464 : Assert( m_ndst < m_mytransfer.size(), "Indexing out of mytransfer[dst]" );
[ - - ][ - - ]
350 [ + + ]: 2464 : if (toMesh == m_meshid) {
351 [ + - ][ + - ]: 1232 : exam2m::setDestPoints( thisProxy, thisIndex, &m_coord, u,
352 : : cb_xfer );
353 : 1232 : ++m_ndst;
354 : : } else {
355 : 1232 : m_ndst = 0;
356 : : }
357 : :
358 : : }
359 : :
360 : 12264 : m_nsrc = 0;
361 : 12264 : m_ndst = 0;
362 : 12264 : }
363 : :
364 : 616 : void Discretization::to_complete()
365 : : // *****************************************************************************
366 : : //! Solution transfer from background to overset mesh completed (from ExaM2M)
367 : : //! \brief This is called by ExaM2M on the destination mesh when the
368 : : //! transfer completes. Since this is called only on the destination, we find
369 : : //! and notify the corresponding source of the completion.
370 : : // *****************************************************************************
371 : : {
372 : : // Lookup the source disc and notify it of completion
373 [ + + ]: 1232 : for (auto& t : m_transfer) {
374 [ + - ]: 616 : if (m_meshid == t.dst) {
375 [ + - ][ + - ]: 616 : m_disc[ t.src ][ thisIndex ].transfer_complete();
376 : : }
377 : : }
378 : :
379 [ + - ]: 616 : thisProxy[ thisIndex ].transfer_complete();
380 : 616 : }
381 : :
382 : 616 : void Discretization::from_complete()
383 : : // *****************************************************************************
384 : : //! Solution transfer from overset to background mesh completed (from ExaM2M)
385 : : //! \brief This is called by ExaM2M on the destination mesh when the
386 : : //! transfer completes. Since this is called only on the destination, we find
387 : : //! and notify the corresponding source of the completion.
388 : : // *****************************************************************************
389 : : {
390 : : // Lookup the source disc and notify it of completion
391 [ + + ]: 1232 : for (auto& t : m_transfer) {
392 [ + - ]: 616 : if (m_meshid == t.src) {
393 [ + - ][ + - ]: 616 : m_disc[ t.dst ][ thisIndex ].transfer_complete_from_dest();
394 : : }
395 : : }
396 : :
397 : 616 : m_transfer_complete.send();
398 : 616 : }
399 : :
400 : 616 : void Discretization::transfer_complete_from_dest()
401 : : // *****************************************************************************
402 : : //! Solution transfer completed (from dest Discretization)
403 : : //! \details Called on the source only by the destination when a back and forth
404 : : //! transfer step completes.
405 : : // *****************************************************************************
406 : : {
407 : 616 : m_transfer_complete.send();
408 : 616 : }
409 : :
410 : 1232 : void Discretization::transfer_complete()
411 : : // *****************************************************************************
412 : : //! Solution transfer completed (one-way)
413 : : //! \note Single exit point after solution transfer between meshes
414 : : // *****************************************************************************
415 : : {
416 [ + - ]: 1232 : contribute( sizeof(nullptr), nullptr, CkReduction::nop,
417 : 2464 : CkCallback(CkReductionTarget(Transporter,solutionTransferred),
418 : : m_transporter) );
419 : 1232 : }
420 : :
421 : : std::vector< std::size_t >
422 : 8329 : Discretization::bndel() const
423 : : // *****************************************************************************
424 : : // Find elements along our mesh chunk boundary
425 : : //! \return List of local element ids that have at least a single node
426 : : //! contributing to a chare boundary
427 : : // *****************************************************************************
428 : : {
429 : : // Lambda to find out if a mesh node is shared with another chare
430 : 10327604 : auto shared = [this]( std::size_t i ){
431 [ + + ]: 29309186 : for (const auto& [c,n] : m_nodeCommMap)
432 [ + - ][ + + ]: 21061788 : if (n.find(i) != end(n)) return true;
433 : 8247398 : return false;
434 : 8329 : };
435 : :
436 : : // Find elements along our mesh chunk boundary
437 : 8329 : std::vector< std::size_t > e;
438 [ + + ]: 10335933 : for (std::size_t n=0; n<m_inpoel.size(); ++n)
439 [ + - ][ + + ]: 10327604 : if (shared( m_gid[ m_inpoel[n] ] )) e.push_back( n/4 );
[ + - ]
440 [ + - ]: 8329 : tk::unique( e );
441 : :
442 : 16658 : return e;
443 : : }
444 : :
445 : : void
446 : 0 : Discretization::resizePostAMR(
447 : : const tk::UnsMesh::Chunk& chunk,
448 : : const tk::UnsMesh::Coords& coord,
449 : : const std::unordered_map< std::size_t, std::size_t >& /*amrNodeMap*/,
450 : : const tk::NodeCommMap& nodeCommMap,
451 : : const std::set< std::size_t >& /*removedNodes*/,
452 : : const std::unordered_map< std::size_t, std::set< std::size_t > >& elemblockid )
453 : : // *****************************************************************************
454 : : // Resize mesh data structures after mesh refinement
455 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
456 : : //! \param[in] coord New mesh node coordinates
457 : : //! \param[in] amrNodeMap Node id map after amr (local ids)
458 : : //! \param[in] nodeCommMap New node communication map
459 : : //! \param[in] removedNodes Newly removed mesh node local ids
460 : : //! \param[in] elemblockid New local tet ids associated with mesh block ids
461 : : // *****************************************************************************
462 : : {
463 : 0 : m_el = chunk; // updates m_inpoel, m_gid, m_lid
464 : 0 : m_nodeCommMap.clear();
465 : 0 : m_nodeCommMap = nodeCommMap; // update node communication map
466 : 0 : m_elemblockid.clear();
467 : 0 : m_elemblockid = elemblockid;
468 : :
469 : : // Update mesh volume container size
470 [ - - ]: 0 : m_vol.resize( m_gid.size(), 0.0 );
471 [ - - ][ - - ]: 0 : if (!m_voln.empty()) m_voln.resize( m_gid.size(), 0.0 );
472 : :
473 : : // Regenerate bid data
474 : 0 : tk::destroy(m_bid);
475 : 0 : m_bid = genBid();
476 : :
477 : : // update mesh node coordinates
478 : 0 : m_coord = coord;
479 : :
480 : : // we are no longer during setup
481 : 0 : m_initial = 0;
482 : 0 : }
483 : :
484 : : void
485 : 3023 : Discretization::startvol()
486 : : // *****************************************************************************
487 : : // Get ready for (re-)computing/communicating nodal volumes
488 : : // *****************************************************************************
489 : : {
490 : 3023 : m_nvol = 0;
491 [ + - ]: 3023 : thisProxy[ thisIndex ].wait4vol();
492 : :
493 : : // Zero out mesh volume container
494 [ + - ]: 3023 : std::fill( begin(m_vol), end(m_vol), 0.0 );
495 : :
496 : : // Clear receive buffer that will be used for collecting nodal volumes
497 : 3023 : m_volc.clear();
498 : 3023 : }
499 : :
500 : : void
501 : 572 : Discretization::registerReducers()
502 : : // *****************************************************************************
503 : : // Configure Charm++ reduction types
504 : : //! \details Since this is a [initnode] routine, see the .ci file, the
505 : : //! Charm++ runtime system executes the routine exactly once on every
506 : : //! logical node early on in the Charm++ init sequence. Must be static as
507 : : //! it is called without an object. See also: Section "Initializations at
508 : : //! Program Startup" at in the Charm++ manual
509 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
510 : : // *****************************************************************************
511 : : {
512 : 572 : PDFMerger = CkReduction::addReducer( tk::mergeUniPDFs );
513 : 572 : }
514 : :
515 : : tk::UnsMesh::Coords
516 : 1943 : Discretization::setCoord( const tk::UnsMesh::CoordMap& coordmap )
517 : : // *****************************************************************************
518 : : // Set mesh coordinates based on coordinates map
519 : : // *****************************************************************************
520 : : {
521 [ - + ][ - - ]: 1943 : Assert( coordmap.size() == m_gid.size(), "Size mismatch" );
[ - - ][ - - ]
522 [ - + ][ - - ]: 1943 : Assert( coordmap.size() == m_lid.size(), "Size mismatch" );
[ - - ][ - - ]
523 : :
524 : 1943 : tk::UnsMesh::Coords coord;
525 [ + - ]: 1943 : coord[0].resize( coordmap.size() );
526 [ + - ]: 1943 : coord[1].resize( coordmap.size() );
527 [ + - ]: 1943 : coord[2].resize( coordmap.size() );
528 : :
529 [ + + ]: 187866 : for (const auto& [ gid, coords ] : coordmap) {
530 [ + - ]: 185923 : auto i = tk::cref_find( m_lid, gid );
531 : 185923 : coord[0][i] = coords[0];
532 : 185923 : coord[1][i] = coords[1];
533 : 185923 : coord[2][i] = coords[2];
534 : : }
535 : :
536 : 1943 : return coord;
537 : : }
538 : :
539 : : void
540 : 10 : Discretization::remap(
541 : : const std::unordered_map< std::size_t, std::size_t >& map )
542 : : // *****************************************************************************
543 : : // Remap mesh data based on new local ids
544 : : //! \param[in] map Mapping of old->new local ids
545 : : // *****************************************************************************
546 : : {
547 : : // Remap connectivity containing local IDs
548 [ + + ][ + - ]: 11914 : for (auto& l : m_inpoel) l = tk::cref_find(map,l);
549 : :
550 : : // Remap global->local id map
551 [ + + ][ + - ]: 1104 : for (auto& [g,l] : m_lid) l = tk::cref_find(map,l);
552 : :
553 : : // Remap global->local id map
554 : 10 : auto maxid = std::numeric_limits< std::size_t >::max();
555 [ + - ]: 20 : std::vector< std::size_t > newgid( m_gid.size(), maxid );
556 [ + + ]: 1104 : for (const auto& [o,n] : map) newgid[n] = m_gid[o];
557 : 10 : m_gid = std::move( newgid );
558 : :
559 [ + - ][ - + ]: 1104 : Assert( std::all_of( m_gid.cbegin(), m_gid.cend(),
[ - - ][ - - ]
[ - - ]
560 : : [=](std::size_t i){ return i < maxid; } ),
561 : : "Not all gid have been remapped" );
562 : :
563 : : // Remap nodal volumes (with contributions along chare-boundaries)
564 [ + - ]: 20 : std::vector< tk::real > newvol( m_vol.size(), 0.0 );
565 [ + + ]: 1104 : for (const auto& [o,n] : map) newvol[n] = m_vol[o];
566 : 10 : m_vol = std::move( newvol );
567 : :
568 : : // Remap nodal volumes (without contributions along chare-boundaries)
569 [ + - ]: 20 : std::vector< tk::real > newv( m_v.size(), 0.0 );
570 [ + + ]: 1104 : for (const auto& [o,n] : map) newv[n] = m_v[o];
571 : 10 : m_v = std::move( newv );
572 : :
573 : : // Remap locations of node coordinates
574 : 20 : tk::UnsMesh::Coords newcoord;
575 : 10 : auto npoin = m_coord[0].size();
576 [ + - ]: 10 : newcoord[0].resize( npoin );
577 [ + - ]: 10 : newcoord[1].resize( npoin );
578 [ + - ]: 10 : newcoord[2].resize( npoin );
579 [ + + ]: 1104 : for (const auto& [o,n] : map) {
580 : 1094 : newcoord[0][n] = m_coord[0][o];
581 : 1094 : newcoord[1][n] = m_coord[1][o];
582 : 1094 : newcoord[2][n] = m_coord[2][o];
583 : : }
584 : 10 : m_coord = std::move( newcoord );
585 : 10 : }
586 : :
587 : : void
588 : 1943 : Discretization::setRefiner( const CProxy_Refiner& ref )
589 : : // *****************************************************************************
590 : : // Set Refiner Charm++ proxy
591 : : //! \param[in] ref Incoming refiner proxy to store
592 : : // *****************************************************************************
593 : : {
594 : 1943 : m_refiner = ref;
595 : 1943 : }
596 : :
597 : : void
598 : 3023 : Discretization::vol()
599 : : // *****************************************************************************
600 : : // Sum mesh volumes to nodes, start communicating them on chare-boundaries
601 : : // *****************************************************************************
602 : : {
603 : 3023 : const auto& x = m_coord[0];
604 : 3023 : const auto& y = m_coord[1];
605 : 3023 : const auto& z = m_coord[2];
606 : :
607 : : // Compute nodal volumes on our chunk of the mesh
608 [ + + ]: 5555430 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
609 : 5552407 : const std::array< std::size_t, 4 > N{{ m_inpoel[e*4+0], m_inpoel[e*4+1],
610 : 5552407 : m_inpoel[e*4+2], m_inpoel[e*4+3] }};
611 : : // compute element Jacobi determinant * 5/120 = element volume / 4
612 : : const std::array< tk::real, 3 >
613 : 5552407 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
614 : 5552407 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
615 : 5552407 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
616 : 5552407 : const auto J = tk::triple( ba, ca, da ) * 5.0 / 120.0;
617 [ - + ][ - - ]: 5552407 : ErrChk( J > 0, "Element Jacobian non-positive: PE:" +
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
618 : : std::to_string(CkMyPe()) + ", node IDs: " +
619 : : std::to_string(m_gid[N[0]]) + ',' +
620 : : std::to_string(m_gid[N[1]]) + ',' +
621 : : std::to_string(m_gid[N[2]]) + ',' +
622 : : std::to_string(m_gid[N[3]]) + ", coords: (" +
623 : : std::to_string(x[N[0]]) + ", " +
624 : : std::to_string(y[N[0]]) + ", " +
625 : : std::to_string(z[N[0]]) + "), (" +
626 : : std::to_string(x[N[1]]) + ", " +
627 : : std::to_string(y[N[1]]) + ", " +
628 : : std::to_string(z[N[1]]) + "), (" +
629 : : std::to_string(x[N[2]]) + ", " +
630 : : std::to_string(y[N[2]]) + ", " +
631 : : std::to_string(z[N[2]]) + "), (" +
632 : : std::to_string(x[N[3]]) + ", " +
633 : : std::to_string(y[N[3]]) + ", " +
634 : : std::to_string(z[N[3]]) + ')' );
635 : : // scatter add V/4 to nodes
636 [ + + ]: 27762035 : for (std::size_t j=0; j<4; ++j) m_vol[N[j]] += J;
637 : :
638 : : // save element volumes at t=t0
639 [ + + ]: 5552407 : if (m_it == 0) m_vol0[e] = J * 4.0;
640 : : }
641 : :
642 : : // Store nodal volumes without contributions from other chares on
643 : : // chare-boundaries
644 : 3023 : m_v = m_vol;
645 : :
646 : : // Send our nodal volume contributions to neighbor chares
647 [ + + ]: 3023 : if (m_nodeCommMap.empty())
648 : 214 : totalvol();
649 : : else
650 [ + + ]: 25341 : for (const auto& [c,n] : m_nodeCommMap) {
651 [ + - ]: 22532 : std::vector< tk::real > v( n.size() );
652 : 22532 : std::size_t j = 0;
653 [ + + ][ + - ]: 437036 : for (auto i : n) v[ j++ ] = m_vol[ tk::cref_find(m_lid,i) ];
654 [ + - ][ + - ]: 22532 : thisProxy[c].comvol( std::vector<std::size_t>(begin(n), end(n)), v );
[ + - ]
655 : : }
656 : :
657 : 3023 : ownvol_complete();
658 : 3023 : }
659 : :
660 : : void
661 : 22532 : Discretization::comvol( const std::vector< std::size_t >& gid,
662 : : const std::vector< tk::real >& nodevol )
663 : : // *****************************************************************************
664 : : // Receive nodal volumes on chare-boundaries
665 : : //! \param[in] gid Global mesh node IDs at which we receive volume contributions
666 : : //! \param[in] nodevol Partial sums of nodal volume contributions to
667 : : //! chare-boundary nodes
668 : : //! \details This function receives contributions to m_vol, which stores the
669 : : //! nodal volumes. While m_vol stores own contributions, m_volc collects the
670 : : //! neighbor chare contributions during communication. This way work on m_vol
671 : : //! and m_volc is overlapped. The contributions are applied in totalvol().
672 : : // *****************************************************************************
673 : : {
674 [ - + ][ - - ]: 22532 : Assert( nodevol.size() == gid.size(), "Size mismatch" );
[ - - ][ - - ]
675 : :
676 [ + + ]: 437036 : for (std::size_t i=0; i<gid.size(); ++i)
677 : 414504 : m_volc[ gid[i] ] += nodevol[i];
678 : :
679 [ + + ]: 22532 : if (++m_nvol == m_nodeCommMap.size()) {
680 : 2809 : m_nvol = 0;
681 : 2809 : comvol_complete();
682 : : }
683 : 22532 : }
684 : :
685 : : void
686 : 3023 : Discretization::totalvol()
687 : : // *****************************************************************************
688 : : // Sum mesh volumes and contribute own mesh volume to total volume
689 : : // *****************************************************************************
690 : : {
691 : : // Add received contributions to nodal volumes
692 [ + + ]: 318860 : for (const auto& [gid, vol] : m_volc)
693 [ + - ]: 315837 : m_vol[ tk::cref_find(m_lid,gid) ] += vol;
694 : :
695 : : // Clear receive buffer
696 : 3023 : tk::destroy(m_volc);
697 : :
698 : : // Sum mesh volume to host
699 : : std::vector< tk::real > tvol{ 0.0,
700 : 3023 : static_cast<tk::real>(m_initial),
701 [ + - ]: 3023 : static_cast<tk::real>(m_meshid) };
702 [ + + ]: 1484616 : for (auto v : m_v) tvol[0] += v;
703 [ + - ]: 3023 : contribute( tvol, CkReduction::sum_double,
704 [ + - ][ + - ]: 6046 : CkCallback(CkReductionTarget(Transporter,totalvol), m_transporter) );
705 : 3023 : }
706 : :
707 : : void
708 : 1943 : Discretization::stat( tk::real mesh_volume )
709 : : // *****************************************************************************
710 : : // Compute mesh cell statistics
711 : : //! \param[in] mesh_volume Total mesh volume
712 : : // *****************************************************************************
713 : : {
714 : : // Store total mesh volume
715 : 1943 : m_meshvol = mesh_volume;
716 : :
717 : 1943 : const auto& x = m_coord[0];
718 : 1943 : const auto& y = m_coord[1];
719 : 1943 : const auto& z = m_coord[2];
720 : :
721 : 1943 : auto MIN = -std::numeric_limits< tk::real >::max();
722 : 1943 : auto MAX = std::numeric_limits< tk::real >::max();
723 [ + - ]: 3886 : std::vector< tk::real > min{ MAX, MAX, MAX };
724 [ + - ]: 3886 : std::vector< tk::real > max{ MIN, MIN, MIN };
725 [ + - ]: 3886 : std::vector< tk::real > sum{ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
726 : 3886 : tk::UniPDF edgePDF( 1e-4 );
727 : 3886 : tk::UniPDF volPDF( 1e-4 );
728 : 3886 : tk::UniPDF ntetPDF( 1e-4 );
729 : :
730 : : // Compute points surrounding points
731 [ + - ][ + - ]: 3886 : auto psup = tk::genPsup( m_inpoel, 4, tk::genEsup(m_inpoel,4) );
732 [ - + ][ - - ]: 1943 : Assert( psup.second.size()-1 == m_gid.size(),
[ - - ][ - - ]
733 : : "Number of mesh points and number of global IDs unequal" );
734 : :
735 : : // Compute edge length statistics
736 : : // Note that while the min and max edge lengths are independent of the number
737 : : // of CPUs (by the time they are aggregated across all chares), the sum of
738 : : // the edge lengths and the edge length PDF are not. This is because the
739 : : // edges on the chare-boundary are counted multiple times and we
740 : : // conscientiously do not make an effort to precisely compute this, because
741 : : // that would require communication and more complex logic. Since these
742 : : // statistics are intended as simple average diagnostics, we ignore these
743 : : // small differences. For reproducible average edge lengths and edge length
744 : : // PDFs, run the mesh in serial.
745 [ + + ]: 187866 : for (std::size_t p=0; p<m_gid.size(); ++p)
746 [ + + ]: 1949919 : for (auto i : tk::Around(psup,p)) {
747 : 1763996 : const auto dx = x[ i ] - x[ p ];
748 : 1763996 : const auto dy = y[ i ] - y[ p ];
749 : 1763996 : const auto dz = z[ i ] - z[ p ];
750 : 1763996 : const auto length = std::sqrt( dx*dx + dy*dy + dz*dz );
751 [ + + ]: 1763996 : if (length < min[0]) min[0] = length;
752 [ + + ]: 1763996 : if (length > max[0]) max[0] = length;
753 : 1763996 : sum[0] += 1.0;
754 : 1763996 : sum[1] += length;
755 [ + - ]: 1763996 : edgePDF.add( length );
756 : : }
757 : :
758 : : // Compute mesh cell volume statistics
759 [ + + ]: 546990 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
760 : 545047 : const std::array< std::size_t, 4 > N{{ m_inpoel[e*4+0], m_inpoel[e*4+1],
761 : 545047 : m_inpoel[e*4+2], m_inpoel[e*4+3] }};
762 : : const std::array< tk::real, 3 >
763 : 545047 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
764 : 545047 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
765 : 545047 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
766 : 545047 : const auto L = std::cbrt( tk::triple( ba, ca, da ) / 6.0 );
767 [ + + ]: 545047 : if (L < min[1]) min[1] = L;
768 [ + + ]: 545047 : if (L > max[1]) max[1] = L;
769 : 545047 : sum[2] += 1.0;
770 : 545047 : sum[3] += L;
771 [ + - ]: 545047 : volPDF.add( L );
772 : : }
773 : :
774 : : // Contribute stats of number of tetrahedra (ntets)
775 : 1943 : sum[4] = 1.0;
776 : 1943 : min[2] = max[2] = sum[5] = static_cast< tk::real >( m_inpoel.size() / 4 );
777 [ + - ]: 1943 : ntetPDF.add( min[2] );
778 : :
779 [ + - ]: 1943 : min.push_back( static_cast<tk::real>(m_meshid) );
780 [ + - ]: 1943 : max.push_back( static_cast<tk::real>(m_meshid) );
781 [ + - ]: 1943 : sum.push_back( static_cast<tk::real>(m_meshid) );
782 : :
783 : : // Contribute to mesh statistics across all Discretization chares
784 [ + - ]: 1943 : contribute( min, CkReduction::min_double,
785 [ + - ][ + - ]: 3886 : CkCallback(CkReductionTarget(Transporter,minstat), m_transporter) );
786 [ + - ]: 1943 : contribute( max, CkReduction::max_double,
787 [ + - ][ + - ]: 3886 : CkCallback(CkReductionTarget(Transporter,maxstat), m_transporter) );
788 [ + - ]: 1943 : contribute( sum, CkReduction::sum_double,
789 [ + - ][ + - ]: 3886 : CkCallback(CkReductionTarget(Transporter,sumstat), m_transporter) );
790 : :
791 : : // Serialize PDFs to raw stream
792 [ + - ][ + - ]: 11658 : auto stream = tk::serialize( m_meshid, { edgePDF, volPDF, ntetPDF } );
[ + - ][ + - ]
[ + - ]
793 : : // Create Charm++ callback function for reduction of PDFs with
794 : : // Transporter::pdfstat() as the final target where the results will appear.
795 [ + - ][ + - ]: 3886 : CkCallback cb( CkIndex_Transporter::pdfstat(nullptr), m_transporter );
796 : : // Contribute serialized PDF of partial sums to host via Charm++ reduction
797 [ + - ]: 1943 : contribute( stream.first, stream.second.get(), PDFMerger, cb );
798 : 1943 : }
799 : :
800 : : void
801 : 1943 : Discretization::boxvol(
802 : : const std::vector< std::unordered_set< std::size_t > >& nodes,
803 : : const std::unordered_map< std::size_t, std::set< std::size_t > >& nodeblk,
804 : : std::size_t nuserblk )
805 : : // *****************************************************************************
806 : : // Compute total box IC volume
807 : : //! \param[in] nodes Node list contributing to box IC volume (for each IC box)
808 : : //! \param[in] nodeblk Node list associated to mesh blocks contributing to block
809 : : //! volumes (for each IC box)
810 : : //! \param[in] nuserblk Number of user IC mesh blocks
811 : : // *****************************************************************************
812 : : {
813 : : // Compute partial box IC volume (just add up all boxes)
814 : 1943 : tk::real boxvol = 0.0;
815 [ + + ][ + + ]: 4485 : for (const auto& b : nodes) for (auto i : b) boxvol += m_v[i];
816 : :
817 : : // Compute partial IC mesh block volume
818 : 3886 : std::vector< tk::real > blockvols;
819 [ - + ]: 1943 : if (nuserblk > 0) {
820 [ - - ]: 0 : blockvols.resize(nuserblk,0.0);
821 [ - - ]: 0 : for (const auto& [blid, ndset] : nodeblk) {
822 : : // The following if-test makes sure we access volumes only of mesh blocks
823 : : // with user-specified ICs
824 [ - - ]: 0 : if (blid < nuserblk) {
825 [ - - ]: 0 : for (const auto& n : ndset) blockvols[blid] += m_v[n];
826 : : }
827 : : }
828 : : }
829 : :
830 : : // Sum up box IC volume across all chares
831 [ + - ]: 1943 : auto meshdata = blockvols;
832 [ + - ]: 1943 : meshdata.push_back(boxvol);
833 [ + - ]: 1943 : meshdata.push_back(static_cast<tk::real>(m_meshid));
834 [ + - ]: 1943 : contribute( meshdata, CkReduction::sum_double,
835 [ + - ][ + - ]: 3886 : CkCallback(CkReductionTarget(Transporter,boxvol), m_transporter) );
836 : 1943 : }
837 : :
838 : : void
839 : 3557 : Discretization::write(
840 : : const std::vector< std::size_t >& inpoel,
841 : : const tk::UnsMesh::Coords& coord,
842 : : const std::map< int, std::vector< std::size_t > >& bface,
843 : : const std::map< int, std::vector< std::size_t > >& bnode,
844 : : const std::vector< std::size_t >& triinpoel,
845 : : const std::vector< std::string>& elemfieldnames,
846 : : const std::vector< std::string>& nodefieldnames,
847 : : const std::vector< std::string>& elemsurfnames,
848 : : const std::vector< std::string>& nodesurfnames,
849 : : const std::vector< std::vector< tk::real > >& elemfields,
850 : : const std::vector< std::vector< tk::real > >& nodefields,
851 : : const std::vector< std::vector< tk::real > >& elemsurfs,
852 : : const std::vector< std::vector< tk::real > >& nodesurfs,
853 : : CkCallback c )
854 : : // *****************************************************************************
855 : : // Output mesh and fields data (solution dump) to file(s)
856 : : //! \param[in] inpoel Mesh connectivity for the mesh chunk to be written
857 : : //! \param[in] coord Node coordinates of the mesh chunk to be written
858 : : //! \param[in] bface Map of boundary-face lists mapped to corresponding side set
859 : : //! ids for this mesh chunk
860 : : //! \param[in] bnode Map of boundary-node lists mapped to corresponding side set
861 : : //! ids for this mesh chunk
862 : : //! \param[in] triinpoel Interconnectivity of points and boundary-face in this
863 : : //! mesh chunk
864 : : //! \param[in] elemfieldnames Names of element fields to be output to file
865 : : //! \param[in] nodefieldnames Names of node fields to be output to file
866 : : //! \param[in] elemsurfnames Names of elemental surface fields to be output to
867 : : //! file
868 : : //! \param[in] nodesurfnames Names of node surface fields to be output to file
869 : : //! \param[in] elemfields Field data in mesh elements to output to file
870 : : //! \param[in] nodefields Field data in mesh nodes to output to file
871 : : //! \param[in] elemsurfs Surface field data in mesh elements to output to file
872 : : //! \param[in] nodesurfs Surface field data in mesh nodes to output to file
873 : : //! \param[in] c Function to continue with after the write
874 : : //! \details Since m_meshwriter is a Charm++ chare group, it never migrates and
875 : : //! an instance is guaranteed on every PE. We index the first PE on every
876 : : //! logical compute node. In Charm++'s non-SMP mode, a node is the same as a
877 : : //! PE, so the index is the same as CkMyPe(). In SMP mode the index is the
878 : : //! first PE on every logical node. In non-SMP mode this yields one or more
879 : : //! output files per PE with zero or non-zero virtualization, respectively. If
880 : : //! there are multiple chares on a PE, the writes are serialized per PE, since
881 : : //! only a single entry method call can be executed at any given time. In SMP
882 : : //! mode, still the same number of files are output (one per chare), but the
883 : : //! output is serialized through the first PE of each compute node. In SMP
884 : : //! mode, channeling multiple files via a single PE on each node is required
885 : : //! by NetCDF and HDF5, as well as ExodusII, since none of these libraries are
886 : : //! thread-safe.
887 : : // *****************************************************************************
888 : : {
889 : : // If the previous iteration refined (or moved) the mesh or this is called
890 : : // before the first time step, we also output the mesh.
891 : 3557 : bool meshoutput = m_itf == 0 ? true : false;
892 : :
893 : 3557 : auto eps = std::numeric_limits< tk::real >::epsilon();
894 : 3557 : bool fieldoutput = false;
895 : :
896 : : // Output field data only if there is no dump at this physical time yet
897 [ + - ]: 3557 : if (std::abs(m_lastDumpTime - m_t) > eps ) {
898 : 3557 : m_lastDumpTime = m_t;
899 : 3557 : ++m_itf;
900 : 3557 : fieldoutput = true;
901 : : }
902 : :
903 : : // set of sidesets where fieldoutput is required
904 : 3557 : std::set< int > outsets;
905 : 3557 : const auto& osv = g_inputdeck.get< tag::field_output, tag::sideset >();
906 [ + - ]: 3557 : outsets.insert(osv.begin(), osv.end());
907 : :
908 : : m_meshwriter[ CkNodeFirst( CkMyNode() ) ].
909 [ + - ][ + - ]: 7114 : write( m_meshid, meshoutput, fieldoutput, m_itr, m_itf, m_t, thisIndex,
910 : 3557 : g_inputdeck.get< tag::cmd, tag::io, tag::output >(),
911 : : inpoel, coord, bface, bnode, triinpoel, elemfieldnames,
912 : : nodefieldnames, elemsurfnames, nodesurfnames, elemfields, nodefields,
913 : : elemsurfs, nodesurfs, outsets, c );
914 : 3557 : }
915 : :
916 : : void
917 : 40348 : Discretization::setdt( tk::real newdt )
918 : : // *****************************************************************************
919 : : // Set time step size
920 : : //! \param[in] newdt Size of the new time step
921 : : // *****************************************************************************
922 : : {
923 : 40348 : m_dtn = m_dt;
924 : 40348 : m_dt = newdt;
925 : :
926 : : // Truncate the size of last time step
927 : 40348 : const auto term = g_inputdeck.get< tag::term >();
928 [ + + ]: 40348 : if (m_t+m_dt > term) m_dt = term - m_t;
929 : 40348 : }
930 : :
931 : : void
932 : 40348 : Discretization::next()
933 : : // *****************************************************************************
934 : : // Prepare for next step
935 : : // *****************************************************************************
936 : : {
937 : : // Update floor of physics time divided by output interval times
938 : 40348 : const auto eps = std::numeric_limits< tk::real >::epsilon();
939 : 40348 : const auto ft = g_inputdeck.get< tag::field_output, tag::time_interval >();
940 [ + - ]: 40348 : if (ft > eps) m_physFieldFloor = std::floor( m_t / ft );
941 : 40348 : const auto ht = g_inputdeck.get< tag::history_output, tag::time_interval >();
942 [ + - ]: 40348 : if (ht > eps) m_physHistFloor = std::floor( m_t / ht );
943 : :
944 : : // Update floors of physics time divided by output interval times for ranges
945 : 40348 : const auto& rf = g_inputdeck.get< tag::field_output, tag::time_range >();
946 [ - + ]: 40348 : if (!rf.empty()) {
947 [ - - ][ - - ]: 0 : if (m_t > rf[0] and m_t < rf[1])
[ - - ]
948 : 0 : m_rangeFieldFloor = std::floor( m_t / rf[2] );
949 : : }
950 : 40348 : const auto& rh = g_inputdeck.get< tag::history_output, tag::time_range >();
951 [ - + ]: 40348 : if (!rh.empty()) {
952 [ - - ][ - - ]: 0 : if (m_t > rh[0] and m_t < rh[1])
[ - - ]
953 : 0 : m_rangeHistFloor = std::floor( m_t / rh[2] );
954 : : }
955 : :
956 : 40348 : ++m_it;
957 : 40348 : m_t += m_dt;
958 : 40348 : }
959 : :
960 : : void
961 : 1948 : Discretization::grindZero()
962 : : // *****************************************************************************
963 : : // Zero grind-time
964 : : // *****************************************************************************
965 : : {
966 : 1948 : m_prevstatus = std::chrono::high_resolution_clock::now();
967 : :
968 [ + + ][ + + ]: 1948 : if (thisIndex == 0 && m_meshid == 0) {
969 : 189 : const auto verbose = g_inputdeck.get< tag::cmd, tag::verbose >();
970 : : const auto& def =
971 : 189 : g_inputdeck_defaults.get< tag::cmd, tag::io, tag::screen >();
972 [ + - ]: 189 : tk::Print print( g_inputdeck.get< tag::cmd >().logname( def, m_nrestart ),
973 : : verbose ? std::cout : std::clog,
974 [ + - ][ + - ]: 378 : std::ios_base::app );
975 [ + - ][ + - ]: 189 : print.diag( "Starting time stepping ..." );
976 : : }
977 : 1948 : }
978 : :
979 : : bool
980 : 38405 : Discretization::restarted( int nrestart )
981 : : // *****************************************************************************
982 : : // Detect if just returned from a checkpoint and if so, zero timers
983 : : //! \param[in] nrestart Number of times restarted
984 : : //! \return True if restart detected
985 : : // *****************************************************************************
986 : : {
987 : : // Detect if just restarted from checkpoint:
988 : : // nrestart == -1 if there was no checkpoint this step
989 : : // d->Nrestart() == nrestart if there was a checkpoint this step
990 : : // if both false, just restarted from a checkpoint
991 [ + + ][ + - ]: 38405 : bool restarted = nrestart != -1 and m_nrestart != nrestart;
992 : :
993 : : // If just restarted from checkpoint
994 [ + + ]: 38405 : if (restarted) {
995 : : // Update number of restarts
996 : 5 : m_nrestart = nrestart;
997 : : // Start timer measuring time stepping wall clock time
998 : 5 : m_timer.zero();
999 : : // Zero grind-timer
1000 : 5 : grindZero();
1001 : : }
1002 : :
1003 : 38405 : return restarted;
1004 : : }
1005 : :
1006 : : std::string
1007 : 88 : Discretization::histfilename( const std::string& id,
1008 : : std::streamsize precision )
1009 : : // *****************************************************************************
1010 : : // Construct history output filename
1011 : : //! \param[in] id History point id
1012 : : //! \param[in] precision Floating point precision to use for output
1013 : : //! \return History file name
1014 : : // *****************************************************************************
1015 : : {
1016 [ + - ]: 176 : auto of = g_inputdeck.get< tag::cmd, tag::io, tag::output >();
1017 [ + - ]: 176 : std::stringstream ss;
1018 : :
1019 : : auto mid =
1020 [ + + ][ + - ]: 252 : m_disc.size() > 1 ? std::string( '.' + std::to_string(m_meshid) ) : "";
[ + - ][ + - ]
[ + + ][ + + ]
[ - - ][ - - ]
1021 [ + - ][ + - ]: 88 : ss << std::setprecision(static_cast<int>(precision)) << of << mid << ".hist." << id;
[ + - ][ + - ]
[ + - ]
1022 : :
1023 [ + - ]: 176 : return ss.str();
1024 : : }
1025 : :
1026 : : void
1027 : 30 : Discretization::histheader( std::vector< std::string >&& names )
1028 : : // *****************************************************************************
1029 : : // Output headers for time history files (one for each point)
1030 : : //! \param[in] names History output variable names
1031 : : // *****************************************************************************
1032 : : {
1033 [ + + ]: 40 : for (const auto& h : m_histdata) {
1034 : 10 : auto prec = g_inputdeck.get< tag::history_output, tag::precision >();
1035 [ + - ]: 10 : tk::DiagWriter hw( histfilename( h.get< tag::id >(), prec ),
1036 : 10 : g_inputdeck.get< tag::history_output, tag::format >(),
1037 [ + - ]: 30 : prec );
1038 [ + - ]: 10 : hw.header( names );
1039 : : }
1040 : 30 : }
1041 : :
1042 : : void
1043 : 254 : Discretization::history( std::vector< std::vector< tk::real > >&& data )
1044 : : // *****************************************************************************
1045 : : // Output time history for a time step
1046 : : //! \param[in] data Time history data for all variables and equations integrated
1047 : : // *****************************************************************************
1048 : : {
1049 [ - + ][ - - ]: 254 : Assert( data.size() == m_histdata.size(), "Size mismatch" );
[ - - ][ - - ]
1050 : :
1051 : 254 : std::size_t i = 0;
1052 [ + + ]: 332 : for (const auto& h : m_histdata) {
1053 : 78 : auto prec = g_inputdeck.get< tag::history_output, tag::precision >();
1054 [ + - ]: 78 : tk::DiagWriter hw( histfilename( h.get< tag::id >(), prec ),
1055 : 78 : g_inputdeck.get< tag::history_output, tag::format >(),
1056 : : prec,
1057 [ + - ]: 156 : std::ios_base::app );
1058 [ + - ]: 78 : hw.diag( m_it, m_t, m_dt, data[i] );
1059 : 78 : ++i;
1060 : : }
1061 : 254 : }
1062 : :
1063 : : bool
1064 : 39573 : Discretization::fielditer() const
1065 : : // *****************************************************************************
1066 : : // Decide if field output iteration count interval is hit
1067 : : //! \return True if field output iteration count interval is hit
1068 : : // *****************************************************************************
1069 : : {
1070 [ + + ]: 39573 : if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;
1071 : :
1072 : 33157 : return m_it % g_inputdeck.get< tag::field_output, tag::interval >() == 0;
1073 : : }
1074 : :
1075 : : bool
1076 : 36028 : Discretization::fieldtime() const
1077 : : // *****************************************************************************
1078 : : // Decide if field output physics time interval is hit
1079 : : //! \return True if field output physics time interval is hit
1080 : : // *****************************************************************************
1081 : : {
1082 [ + + ]: 36028 : if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;
1083 : :
1084 : 29612 : const auto eps = std::numeric_limits< tk::real >::epsilon();
1085 : 29612 : const auto ft = g_inputdeck.get< tag::field_output, tag::time_interval >();
1086 : :
1087 [ - + ]: 29612 : if (ft < eps) return false;
1088 : :
1089 : 29612 : return std::floor(m_t/ft) - m_physFieldFloor > eps;
1090 : : }
1091 : :
1092 : : bool
1093 : 36006 : Discretization::fieldrange() const
1094 : : // *****************************************************************************
1095 : : // Decide if physics time falls into a field output time range
1096 : : //! \return True if physics time falls into a field output time range
1097 : : // *****************************************************************************
1098 : : {
1099 [ + + ]: 36006 : if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;
1100 : :
1101 : 29590 : const auto eps = std::numeric_limits< tk::real >::epsilon();
1102 : :
1103 : 29590 : bool output = false;
1104 : :
1105 : 29590 : const auto& rf = g_inputdeck.get< tag::field_output, tag::time_range >();
1106 [ - + ]: 29590 : if (!rf.empty()) {
1107 [ - - ][ - - ]: 0 : if (m_t > rf[0] and m_t < rf[1])
[ - - ]
1108 : 0 : output |= std::floor(m_t/rf[2]) - m_rangeFieldFloor > eps;
1109 : : }
1110 : :
1111 : 29590 : return output;
1112 : : }
1113 : :
1114 : : bool
1115 : 42562 : Discretization::histiter() const
1116 : : // *****************************************************************************
1117 : : // Decide if history output iteration count interval is hit
1118 : : //! \return True if history output iteration count interval is hit
1119 : : // *****************************************************************************
1120 : : {
1121 : 42562 : const auto hist = g_inputdeck.get< tag::history_output, tag::interval >();
1122 : 42562 : const auto& hist_points = g_inputdeck.get< tag::history_output, tag::point >();
1123 : :
1124 [ + + ][ + - ]: 42562 : return m_it % hist == 0 and not hist_points.empty();
1125 : : }
1126 : :
1127 : : bool
1128 : 42275 : Discretization::histtime() const
1129 : : // *****************************************************************************
1130 : : // Decide if history output physics time interval is hit
1131 : : //! \return True if history output physics time interval is hit
1132 : : // *****************************************************************************
1133 : : {
1134 [ + + ]: 42275 : if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;
1135 : :
1136 : 32355 : const auto eps = std::numeric_limits< tk::real >::epsilon();
1137 : 32355 : const auto ht = g_inputdeck.get< tag::history_output, tag::time_interval >();
1138 : :
1139 [ - + ]: 32355 : if (ht < eps) return false;
1140 : :
1141 : 32355 : return std::floor(m_t/ht) - m_physHistFloor > eps;
1142 : : }
1143 : :
1144 : : bool
1145 : 42269 : Discretization::histrange() const
1146 : : // *****************************************************************************
1147 : : // Decide if physics time falls into a history output time range
1148 : : //! \return True if physics time falls into a history output time range
1149 : : // *****************************************************************************
1150 : : {
1151 [ + + ]: 42269 : if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;
1152 : :
1153 : 32349 : const auto eps = std::numeric_limits< tk::real >::epsilon();
1154 : :
1155 : 32349 : bool output = false;
1156 : :
1157 : 32349 : const auto& rh = g_inputdeck.get< tag::history_output, tag::time_range >();
1158 [ - + ]: 32349 : if (!rh.empty()) {
1159 [ - - ][ - - ]: 0 : if (m_t > rh[0] and m_t < rh[1])
[ - - ]
1160 : 0 : output |= std::floor(m_t/rh[2]) - m_rangeHistFloor > eps;
1161 : : }
1162 : :
1163 : 32349 : return output;
1164 : : }
1165 : :
1166 : : bool
1167 : 50124 : Discretization::finished() const
1168 : : // *****************************************************************************
1169 : : // Decide if this is the last time step
1170 : : //! \return True if this is the last time step
1171 : : // *****************************************************************************
1172 : : {
1173 : 50124 : const auto eps = std::numeric_limits< tk::real >::epsilon();
1174 : 50124 : const auto nstep = g_inputdeck.get< tag::nstep >();
1175 : 50124 : const auto term = g_inputdeck.get< tag::term >();
1176 : :
1177 [ + + ][ + + ]: 50124 : return std::abs(m_t-term) < eps or m_it >= nstep;
1178 : : }
1179 : :
1180 : : void
1181 : 40348 : Discretization::status()
1182 : : // *****************************************************************************
1183 : : // Output one-liner status report
1184 : : // *****************************************************************************
1185 : : {
1186 : : // Query after how many time steps user wants TTY dump
1187 : 40348 : const auto tty = g_inputdeck.get< tag::ttyi >();
1188 : :
1189 : : // estimate grind time (taken between this and the previous time step)
1190 : : using std::chrono::duration_cast;
1191 : : using ms = std::chrono::milliseconds;
1192 : : using clock = std::chrono::high_resolution_clock;
1193 [ + - ][ + - ]: 40348 : auto grind_time = duration_cast< ms >(clock::now() - m_prevstatus).count();
1194 : 40348 : m_prevstatus = clock::now();
1195 : :
1196 [ + + ][ + + ]: 40348 : if (thisIndex==0 and m_meshid == 0 and not (m_it%tty)) {
[ + + ]
1197 : :
1198 : 2214 : const auto term = g_inputdeck.get< tag::term >();
1199 : 2214 : const auto t0 = g_inputdeck.get< tag::t0 >();
1200 : 2214 : const auto nstep = g_inputdeck.get< tag::nstep >();
1201 : 2214 : const auto diag = g_inputdeck.get< tag::diagnostics, tag::interval >();
1202 : 2214 : const auto lbfreq = g_inputdeck.get< tag::cmd, tag::lbfreq >();
1203 : 2214 : const auto rsfreq = g_inputdeck.get< tag::cmd, tag::rsfreq >();
1204 : 2214 : const auto verbose = g_inputdeck.get< tag::cmd, tag::verbose >();
1205 : 2214 : const auto benchmark = g_inputdeck.get< tag::cmd, tag::benchmark >();
1206 : 2214 : const auto steady = g_inputdeck.get< tag::steady_state >();
1207 : :
1208 : : // estimate time elapsed and time for accomplishment
1209 [ + - ][ + - ]: 2214 : tk::Timer::Watch ete, eta;
1210 [ + + ][ + - ]: 2214 : if (not steady) m_timer.eta( term-t0, m_t-t0, nstep, m_it, ete, eta );
1211 : :
1212 : : const auto& def =
1213 : 2214 : g_inputdeck_defaults.get< tag::cmd, tag::io, tag::screen >();
1214 [ + - ]: 2214 : tk::Print print( g_inputdeck.get< tag::cmd >().logname( def, m_nrestart ),
1215 : : verbose ? std::cout : std::clog,
1216 [ + - ][ + - ]: 6642 : std::ios_base::app );
1217 : :
1218 : : // Output one-liner
1219 : 2214 : print << std::setfill(' ') << std::setw(8) << m_it << " "
1220 : 0 : << std::scientific << std::setprecision(6)
1221 : 0 : << std::setw(12) << m_t << " "
1222 : 2214 : << m_dt << " "
1223 : 0 : << std::setfill('0')
1224 : 0 : << std::setw(3) << ete.hrs.count() << ":"
1225 : 0 : << std::setw(2) << ete.min.count() << ":"
1226 : 0 : << std::setw(2) << ete.sec.count() << " "
1227 : 0 : << std::setw(3) << eta.hrs.count() << ":"
1228 : 0 : << std::setw(2) << eta.min.count() << ":"
1229 : 0 : << std::setw(2) << eta.sec.count() << " "
1230 : 0 : << std::scientific << std::setprecision(6) << std::setfill(' ')
1231 [ + - ][ + - ]: 2214 : << std::setw(9) << grind_time << " ";
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
1232 : :
1233 : : // Augment one-liner status with output indicators
1234 [ + - ][ + + ]: 2214 : if (fielditer() or fieldtime() or fieldrange()) print << 'f';
[ + - ][ + + ]
[ + - ][ - + ]
[ + + ][ + - ]
1235 [ + + ][ + - ]: 2214 : if (not (m_it % diag)) print << 'd';
1236 [ + - ][ + + ]: 2214 : if (histiter() or histtime() or histrange()) print << 't';
[ + - ][ + + ]
[ + - ][ - + ]
[ + + ][ + - ]
1237 [ - + ][ - - ]: 2214 : if (m_refined) print << 'h';
1238 [ + + ][ + - ]: 2214 : if (not (m_it % lbfreq) && not finished()) print << 'l';
[ + + ][ + + ]
[ + - ]
1239 [ + + ][ + - ]: 2214 : if (not benchmark && (not (m_it % rsfreq) || finished())) print << 'r';
[ + - ][ + + ]
[ + + ][ + - ]
1240 : :
1241 [ + + ][ + - ]: 2214 : if (not m_meshvel_converged) print << 'a';
1242 : 2214 : m_meshvel_converged = true; // get ready for next time step
1243 : :
1244 [ + - ]: 2214 : print << std::endl;
1245 : : }
1246 : 40348 : }
1247 : :
1248 : : #include "NoWarning/discretization.def.h"
|
