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