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 : 1994 : 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 : 1994 : int nc ) :
54 : : m_meshid( meshid ),
55 : 1994 : 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 : 3988 : m_t( g_inputdeck.get< tag::t0 >() ),
63 : 3988 : 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 : 3988 : m_dt( g_inputdeck.get< tag::dt >() ),
69 : 1994 : 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 : 1994 : 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 : 1994 : m_voln( m_vol ),
85 : 1994 : m_vol0( m_inpoel.size()/4, 0.0 ),
86 : : m_bid(),
87 : : m_timer(),
88 : : m_refined( 0 ),
89 : 1994 : 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 [ + - ][ + - ]: 3988 : 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 [ - + ][ - - ]: 1994 : Assert( !m_inpoel.empty(), "No elements assigned to Discretization chare" );
[ - - ][ - - ]
119 [ + - ][ - + ]: 1994 : 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 [ + - ][ - + ]: 1994 : Assert( tk::conforming( m_inpoel, m_coord ),
[ - - ][ - - ]
[ - - ]
129 : : "Input mesh to Discretization not conforming" );
130 : : #endif
131 : :
132 : : // Store communication maps
133 [ + + ]: 21886 : for (const auto& [ c, maps ] : msum) {
134 [ + - ][ + - ]: 19892 : m_nodeCommMap[c] = maps.get< tag::node >();
135 [ + - ][ + - ]: 19892 : m_edgeCommMap[c] = maps.get< tag::edge >();
136 : : }
137 : :
138 : : // Get ready for computing/communicating nodal volumes
139 [ + - ]: 1994 : startvol();
140 : :
141 : : // Get chare-boundary node-id map
142 [ + - ]: 1994 : 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 : 1994 : const auto& pt = g_inputdeck.get< tag::history_output, tag::point >();
147 [ + + ]: 2046 : 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 [ + + ]: 24180 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
152 [ + - ][ + + ]: 24138 : 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 [ + + ]: 1994 : 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 [ + - ]: 1963 : m_meshvel.resize( m_gid.size() );
166 : : }
167 : :
168 : : // Register mesh with mesh-transfer lib
169 [ + + ][ - + ]: 1994 : if (m_disc.size() == 1 || m_transfer.empty()) {
[ + + ]
170 : : // skip transfer if single mesh or if not involved in coupling
171 [ + - ]: 1980 : 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 : 1994 : }
180 : :
181 : : std::unordered_map< std::size_t, std::size_t >
182 : 1994 : 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 [ + - ]: 3988 : std::vector< std::size_t > c( tk::sumvalsize( m_nodeCommMap ) );
190 : 1994 : std::size_t j = 0;
191 [ + + ][ + + ]: 150736 : for (const auto& [ch,n] : m_nodeCommMap) for (auto i : n) c[j++] = i;
192 [ + - ]: 1994 : tk::unique( c );
193 [ + - ]: 3988 : return tk::assignLid( c );
194 : : }
195 : :
196 : : void
197 : 1994 : 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 : 1994 : std::size_t npoin = m_gid.size();
205 [ + + ][ + - ]: 193778 : 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 [ + - ]: 1994 : std::vector< std::size_t > meshdata{ m_meshid, npoin };
210 [ + - ]: 1994 : contribute( meshdata, CkReduction::sum_ulong,
211 [ + - ][ + - ]: 3988 : CkCallback( CkReductionTarget(Transporter,disccreated), m_transporter ) );
212 : 1994 : }
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 [ + + ]: 29293769 : for (const auto& [c,n] : m_nodeCommMap)
432 [ + - ][ + + ]: 21046371 : 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 : 3074 : Discretization::startvol()
486 : : // *****************************************************************************
487 : : // Get ready for (re-)computing/communicating nodal volumes
488 : : // *****************************************************************************
489 : : {
490 : 3074 : m_nvol = 0;
491 [ + - ]: 3074 : thisProxy[ thisIndex ].wait4vol();
492 : :
493 : : // Zero out mesh volume container
494 [ + - ]: 3074 : std::fill( begin(m_vol), end(m_vol), 0.0 );
495 : :
496 : : // Clear receive buffer that will be used for collecting nodal volumes
497 : 3074 : m_volc.clear();
498 : 3074 : }
499 : :
500 : : void
501 : 587 : 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 : 587 : PDFMerger = CkReduction::addReducer( tk::mergeUniPDFs );
513 : 587 : }
514 : :
515 : : tk::UnsMesh::Coords
516 : 1994 : Discretization::setCoord( const tk::UnsMesh::CoordMap& coordmap )
517 : : // *****************************************************************************
518 : : // Set mesh coordinates based on coordinates map
519 : : // *****************************************************************************
520 : : {
521 [ - + ][ - - ]: 1994 : Assert( coordmap.size() == m_gid.size(), "Size mismatch" );
[ - - ][ - - ]
522 [ - + ][ - - ]: 1994 : Assert( coordmap.size() == m_lid.size(), "Size mismatch" );
[ - - ][ - - ]
523 : :
524 : 1994 : tk::UnsMesh::Coords coord;
525 [ + - ]: 1994 : coord[0].resize( coordmap.size() );
526 [ + - ]: 1994 : coord[1].resize( coordmap.size() );
527 [ + - ]: 1994 : coord[2].resize( coordmap.size() );
528 : :
529 [ + + ]: 193778 : for (const auto& [ gid, coords ] : coordmap) {
530 [ + - ]: 191784 : auto i = tk::cref_find( m_lid, gid );
531 : 191784 : coord[0][i] = coords[0];
532 : 191784 : coord[1][i] = coords[1];
533 : 191784 : coord[2][i] = coords[2];
534 : : }
535 : :
536 : 1994 : 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 : 1994 : Discretization::setRefiner( const CProxy_Refiner& ref )
589 : : // *****************************************************************************
590 : : // Set Refiner Charm++ proxy
591 : : //! \param[in] ref Incoming refiner proxy to store
592 : : // *****************************************************************************
593 : : {
594 : 1994 : m_refiner = ref;
595 : 1994 : }
596 : :
597 : : void
598 : 3074 : Discretization::vol()
599 : : // *****************************************************************************
600 : : // Sum mesh volumes to nodes, start communicating them on chare-boundaries
601 : : // *****************************************************************************
602 : : {
603 : 3074 : const auto& x = m_coord[0];
604 : 3074 : const auto& y = m_coord[1];
605 : 3074 : const auto& z = m_coord[2];
606 : :
607 : : // Compute nodal volumes on our chunk of the mesh
608 [ + + ]: 5570338 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
609 : 5567264 : const std::array< std::size_t, 4 > N{{ m_inpoel[e*4+0], m_inpoel[e*4+1],
610 : 5567264 : 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 : 5567264 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
614 : 5567264 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
615 : 5567264 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
616 : 5567264 : const auto J = tk::triple( ba, ca, da ) * 5.0 / 120.0;
617 [ - + ][ - - ]: 5567264 : 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 [ + + ]: 27836320 : for (std::size_t j=0; j<4; ++j) m_vol[N[j]] += J;
637 : :
638 : : // save element volumes at t=t0
639 [ + + ]: 5567264 : 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 : 3074 : m_v = m_vol;
645 : :
646 : : // Send our nodal volume contributions to neighbor chares
647 [ + + ]: 3074 : if (m_nodeCommMap.empty())
648 : 217 : totalvol();
649 : : else
650 [ + + ]: 25629 : for (const auto& [c,n] : m_nodeCommMap) {
651 [ + - ]: 22772 : std::vector< tk::real > v( n.size() );
652 : 22772 : std::size_t j = 0;
653 [ + + ][ + - ]: 439502 : for (auto i : n) v[ j++ ] = m_vol[ tk::cref_find(m_lid,i) ];
654 [ + - ][ + - ]: 22772 : thisProxy[c].comvol( std::vector<std::size_t>(begin(n), end(n)), v );
[ + - ]
655 : : }
656 : :
657 : 3074 : ownvol_complete();
658 : 3074 : }
659 : :
660 : : void
661 : 22772 : 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 [ - + ][ - - ]: 22772 : Assert( nodevol.size() == gid.size(), "Size mismatch" );
[ - - ][ - - ]
675 : :
676 [ + + ]: 439502 : for (std::size_t i=0; i<gid.size(); ++i)
677 : 416730 : m_volc[ gid[i] ] += nodevol[i];
678 : :
679 [ + + ]: 22772 : if (++m_nvol == m_nodeCommMap.size()) {
680 : 2857 : m_nvol = 0;
681 : 2857 : comvol_complete();
682 : : }
683 : 22772 : }
684 : :
685 : : void
686 : 3074 : 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 [ + + ]: 320676 : for (const auto& [gid, vol] : m_volc)
693 [ + - ]: 317602 : m_vol[ tk::cref_find(m_lid,gid) ] += vol;
694 : :
695 : : // Clear receive buffer
696 : 3074 : tk::destroy(m_volc);
697 : :
698 : : // Sum mesh volume to host
699 : : std::vector< tk::real > tvol{ 0.0,
700 : 3074 : static_cast<tk::real>(m_initial),
701 [ + - ]: 3074 : static_cast<tk::real>(m_meshid) };
702 [ + + ]: 1490528 : for (auto v : m_v) tvol[0] += v;
703 [ + - ]: 3074 : contribute( tvol, CkReduction::sum_double,
704 [ + - ][ + - ]: 6148 : CkCallback(CkReductionTarget(Transporter,totalvol), m_transporter) );
705 : 3074 : }
706 : :
707 : : void
708 : 1994 : 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 : 1994 : m_meshvol = mesh_volume;
716 : :
717 : 1994 : const auto& x = m_coord[0];
718 : 1994 : const auto& y = m_coord[1];
719 : 1994 : const auto& z = m_coord[2];
720 : :
721 : 1994 : auto MIN = -std::numeric_limits< tk::real >::max();
722 : 1994 : auto MAX = std::numeric_limits< tk::real >::max();
723 [ + - ]: 3988 : std::vector< tk::real > min{ MAX, MAX, MAX };
724 [ + - ]: 3988 : std::vector< tk::real > max{ MIN, MIN, MIN };
725 [ + - ]: 3988 : std::vector< tk::real > sum{ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
726 : 3988 : tk::UniPDF edgePDF( 1e-4 );
727 : 3988 : tk::UniPDF volPDF( 1e-4 );
728 : 3988 : tk::UniPDF ntetPDF( 1e-4 );
729 : :
730 : : // Compute points surrounding points
731 [ + - ][ + - ]: 3988 : auto psup = tk::genPsup( m_inpoel, 4, tk::genEsup(m_inpoel,4) );
732 [ - + ][ - - ]: 1994 : 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 [ + + ]: 193778 : for (std::size_t p=0; p<m_gid.size(); ++p)
746 [ + + ]: 2008234 : for (auto i : tk::Around(psup,p)) {
747 : 1816450 : const auto dx = x[ i ] - x[ p ];
748 : 1816450 : const auto dy = y[ i ] - y[ p ];
749 : 1816450 : const auto dz = z[ i ] - z[ p ];
750 : 1816450 : const auto length = std::sqrt( dx*dx + dy*dy + dz*dz );
751 [ + + ]: 1816450 : if (length < min[0]) min[0] = length;
752 [ + + ]: 1816450 : if (length > max[0]) max[0] = length;
753 : 1816450 : sum[0] += 1.0;
754 : 1816450 : sum[1] += length;
755 [ + - ]: 1816450 : edgePDF.add( length );
756 : : }
757 : :
758 : : // Compute mesh cell volume statistics
759 [ + + ]: 561898 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
760 : 559904 : const std::array< std::size_t, 4 > N{{ m_inpoel[e*4+0], m_inpoel[e*4+1],
761 : 559904 : m_inpoel[e*4+2], m_inpoel[e*4+3] }};
762 : : const std::array< tk::real, 3 >
763 : 559904 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
764 : 559904 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
765 : 559904 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
766 : 559904 : const auto L = std::cbrt( tk::triple( ba, ca, da ) / 6.0 );
767 [ + + ]: 559904 : if (L < min[1]) min[1] = L;
768 [ + + ]: 559904 : if (L > max[1]) max[1] = L;
769 : 559904 : sum[2] += 1.0;
770 : 559904 : sum[3] += L;
771 [ + - ]: 559904 : volPDF.add( L );
772 : : }
773 : :
774 : : // Contribute stats of number of tetrahedra (ntets)
775 : 1994 : sum[4] = 1.0;
776 : 1994 : min[2] = max[2] = sum[5] = static_cast< tk::real >( m_inpoel.size() / 4 );
777 [ + - ]: 1994 : ntetPDF.add( min[2] );
778 : :
779 [ + - ]: 1994 : min.push_back( static_cast<tk::real>(m_meshid) );
780 [ + - ]: 1994 : max.push_back( static_cast<tk::real>(m_meshid) );
781 [ + - ]: 1994 : sum.push_back( static_cast<tk::real>(m_meshid) );
782 : :
783 : : // Contribute to mesh statistics across all Discretization chares
784 [ + - ]: 1994 : contribute( min, CkReduction::min_double,
785 [ + - ][ + - ]: 3988 : CkCallback(CkReductionTarget(Transporter,minstat), m_transporter) );
786 [ + - ]: 1994 : contribute( max, CkReduction::max_double,
787 [ + - ][ + - ]: 3988 : CkCallback(CkReductionTarget(Transporter,maxstat), m_transporter) );
788 [ + - ]: 1994 : contribute( sum, CkReduction::sum_double,
789 [ + - ][ + - ]: 3988 : CkCallback(CkReductionTarget(Transporter,sumstat), m_transporter) );
790 : :
791 : : // Serialize PDFs to raw stream
792 [ + - ][ + - ]: 11964 : 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 [ + - ][ + - ]: 3988 : CkCallback cb( CkIndex_Transporter::pdfstat(nullptr), m_transporter );
796 : : // Contribute serialized PDF of partial sums to host via Charm++ reduction
797 [ + - ]: 1994 : contribute( stream.first, stream.second.get(), PDFMerger, cb );
798 : 1994 : }
799 : :
800 : : void
801 : 1994 : 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 : 1994 : tk::real boxvol = 0.0;
815 [ + + ][ + + ]: 4536 : for (const auto& b : nodes) for (auto i : b) boxvol += m_v[i];
816 : :
817 : : // Compute partial IC mesh block volume
818 : 3988 : std::vector< tk::real > blockvols;
819 [ - + ]: 1994 : 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 [ + - ]: 1994 : auto meshdata = blockvols;
832 [ + - ]: 1994 : meshdata.push_back(boxvol);
833 [ + - ]: 1994 : meshdata.push_back(static_cast<tk::real>(m_meshid));
834 [ + - ]: 1994 : contribute( meshdata, CkReduction::sum_double,
835 [ + - ][ + - ]: 3988 : CkCallback(CkReductionTarget(Transporter,boxvol), m_transporter) );
836 : 1994 : }
837 : :
838 : : void
839 : 3709 : 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 : 3709 : bool meshoutput = m_itf == 0 ? true : false;
892 : :
893 : 3709 : auto eps = std::numeric_limits< tk::real >::epsilon();
894 : 3709 : bool fieldoutput = false;
895 : :
896 : : // Output field data only if there is no dump at this physical time yet
897 [ + - ]: 3709 : if (std::abs(m_lastDumpTime - m_t) > eps ) {
898 : 3709 : m_lastDumpTime = m_t;
899 : 3709 : ++m_itf;
900 : 3709 : fieldoutput = true;
901 : : }
902 : :
903 : : // set of sidesets where fieldoutput is required
904 : 3709 : std::set< int > outsets;
905 : 3709 : const auto& osv = g_inputdeck.get< tag::field_output, tag::sideset >();
906 [ + - ]: 3709 : outsets.insert(osv.begin(), osv.end());
907 : :
908 : : m_meshwriter[ CkNodeFirst( CkMyNode() ) ].
909 [ + - ][ + - ]: 7418 : write( m_meshid, meshoutput, fieldoutput, m_itr, m_itf, m_t, thisIndex,
910 : 3709 : 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 : 3709 : }
915 : :
916 : : void
917 : 42748 : Discretization::setdt( tk::real newdt )
918 : : // *****************************************************************************
919 : : // Set time step size
920 : : //! \param[in] newdt Size of the new time step
921 : : // *****************************************************************************
922 : : {
923 : 42748 : m_dtn = m_dt;
924 : 42748 : m_dt = newdt;
925 : :
926 : : // Truncate the size of last time step
927 : 42748 : const auto term = g_inputdeck.get< tag::term >();
928 [ + + ]: 42748 : if (m_t+m_dt > term) m_dt = term - m_t;
929 : 42748 : }
930 : :
931 : : void
932 : 42748 : Discretization::next()
933 : : // *****************************************************************************
934 : : // Prepare for next step
935 : : // *****************************************************************************
936 : : {
937 : : // Update floor of physics time divided by output interval times
938 : 42748 : const auto eps = std::numeric_limits< tk::real >::epsilon();
939 : 42748 : const auto ft = g_inputdeck.get< tag::field_output, tag::time_interval >();
940 [ + - ]: 42748 : if (ft > eps) m_physFieldFloor = std::floor( m_t / ft );
941 : 42748 : const auto ht = g_inputdeck.get< tag::history_output, tag::time_interval >();
942 [ + - ]: 42748 : 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 : 42748 : const auto& rf = g_inputdeck.get< tag::field_output, tag::time_range >();
946 [ - + ]: 42748 : 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 : 0 : const auto& rh = g_inputdeck.get< tag::history_output, tag::time_range >();
950 [ - - ][ - - ]: 0 : if (m_t > rh[0] and m_t < rh[1])
[ - - ]
951 : 0 : m_rangeHistFloor = std::floor( m_t / rh[2] );
952 : : }
953 : :
954 : 42748 : ++m_it;
955 : 42748 : m_t += m_dt;
956 : 42748 : }
957 : :
958 : : void
959 : 1999 : Discretization::grindZero()
960 : : // *****************************************************************************
961 : : // Zero grind-time
962 : : // *****************************************************************************
963 : : {
964 : 1999 : m_prevstatus = std::chrono::high_resolution_clock::now();
965 : :
966 [ + + ][ + + ]: 1999 : if (thisIndex == 0 && m_meshid == 0) {
967 : 195 : const auto verbose = g_inputdeck.get< tag::cmd, tag::verbose >();
968 : : const auto& def =
969 : 195 : g_inputdeck_defaults.get< tag::cmd, tag::io, tag::screen >();
970 [ + - ]: 195 : tk::Print print( g_inputdeck.get< tag::cmd >().logname( def, m_nrestart ),
971 : : verbose ? std::cout : std::clog,
972 [ + - ][ + - ]: 390 : std::ios_base::app );
973 [ + - ][ + - ]: 195 : print.diag( "Starting time stepping ..." );
974 : : }
975 : 1999 : }
976 : :
977 : : bool
978 : 40754 : Discretization::restarted( int nrestart )
979 : : // *****************************************************************************
980 : : // Detect if just returned from a checkpoint and if so, zero timers
981 : : //! \param[in] nrestart Number of times restarted
982 : : //! \return True if restart detected
983 : : // *****************************************************************************
984 : : {
985 : : // Detect if just restarted from checkpoint:
986 : : // nrestart == -1 if there was no checkpoint this step
987 : : // d->Nrestart() == nrestart if there was a checkpoint this step
988 : : // if both false, just restarted from a checkpoint
989 [ + + ][ + - ]: 40754 : bool restarted = nrestart != -1 and m_nrestart != nrestart;
990 : :
991 : : // If just restarted from checkpoint
992 [ + + ]: 40754 : if (restarted) {
993 : : // Update number of restarts
994 : 5 : m_nrestart = nrestart;
995 : : // Start timer measuring time stepping wall clock time
996 : 5 : m_timer.zero();
997 : : // Zero grind-timer
998 : 5 : grindZero();
999 : : }
1000 : :
1001 : 40754 : return restarted;
1002 : : }
1003 : :
1004 : : std::string
1005 : 88 : Discretization::histfilename( const std::string& id,
1006 : : std::streamsize precision )
1007 : : // *****************************************************************************
1008 : : // Construct history output filename
1009 : : //! \param[in] id History point id
1010 : : //! \param[in] precision Floating point precision to use for output
1011 : : //! \return History file name
1012 : : // *****************************************************************************
1013 : : {
1014 [ + - ]: 176 : auto of = g_inputdeck.get< tag::cmd, tag::io, tag::output >();
1015 [ + - ]: 176 : std::stringstream ss;
1016 : :
1017 : : auto mid =
1018 [ + + ][ + - ]: 252 : m_disc.size() > 1 ? std::string( '.' + std::to_string(m_meshid) ) : "";
[ + - ][ + - ]
[ + + ][ + + ]
[ - - ][ - - ]
1019 [ + - ][ + - ]: 88 : ss << std::setprecision(static_cast<int>(precision)) << of << mid << ".hist." << id;
[ + - ][ + - ]
[ + - ]
1020 : :
1021 [ + - ]: 176 : return ss.str();
1022 : : }
1023 : :
1024 : : void
1025 : 30 : Discretization::histheader( std::vector< std::string >&& names )
1026 : : // *****************************************************************************
1027 : : // Output headers for time history files (one for each point)
1028 : : //! \param[in] names History output variable names
1029 : : // *****************************************************************************
1030 : : {
1031 [ + + ]: 40 : for (const auto& h : m_histdata) {
1032 : 10 : auto prec = g_inputdeck.get< tag::history_output, tag::precision >();
1033 [ + - ]: 10 : tk::DiagWriter hw( histfilename( h.get< tag::id >(), prec ),
1034 : 10 : g_inputdeck.get< tag::history_output, tag::format >(),
1035 [ + - ]: 30 : prec );
1036 [ + - ]: 10 : hw.header( names );
1037 : : }
1038 : 30 : }
1039 : :
1040 : : void
1041 : 254 : Discretization::history( std::vector< std::vector< tk::real > >&& data )
1042 : : // *****************************************************************************
1043 : : // Output time history for a time step
1044 : : //! \param[in] data Time history data for all variables and equations integrated
1045 : : // *****************************************************************************
1046 : : {
1047 [ - + ][ - - ]: 254 : Assert( data.size() == m_histdata.size(), "Size mismatch" );
[ - - ][ - - ]
1048 : :
1049 : 254 : std::size_t i = 0;
1050 [ + + ]: 332 : for (const auto& h : m_histdata) {
1051 : 78 : auto prec = g_inputdeck.get< tag::history_output, tag::precision >();
1052 [ + - ]: 78 : tk::DiagWriter hw( histfilename( h.get< tag::id >(), prec ),
1053 : 78 : g_inputdeck.get< tag::history_output, tag::format >(),
1054 : : prec,
1055 [ + - ]: 156 : std::ios_base::app );
1056 [ + - ]: 78 : hw.diag( m_it, m_t, m_dt, data[i] );
1057 : 78 : ++i;
1058 : : }
1059 : 254 : }
1060 : :
1061 : : bool
1062 : 42196 : Discretization::fielditer() const
1063 : : // *****************************************************************************
1064 : : // Decide if field output iteration count interval is hit
1065 : : //! \return True if field output iteration count interval is hit
1066 : : // *****************************************************************************
1067 : : {
1068 [ + + ]: 42196 : if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;
1069 : :
1070 : 35780 : return m_it % g_inputdeck.get< tag::field_output, tag::interval >() == 0;
1071 : : }
1072 : :
1073 : : bool
1074 : 38488 : Discretization::fieldtime() const
1075 : : // *****************************************************************************
1076 : : // Decide if field output physics time interval is hit
1077 : : //! \return True if field output physics time interval is hit
1078 : : // *****************************************************************************
1079 : : {
1080 [ + + ]: 38488 : if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;
1081 : :
1082 : 32072 : const auto eps = std::numeric_limits< tk::real >::epsilon();
1083 : 32072 : const auto ft = g_inputdeck.get< tag::field_output, tag::time_interval >();
1084 : :
1085 [ - + ]: 32072 : if (ft < eps) return false;
1086 : :
1087 : 32072 : return std::floor(m_t/ft) - m_physFieldFloor > eps;
1088 : : }
1089 : :
1090 : : bool
1091 : 38466 : Discretization::fieldrange() const
1092 : : // *****************************************************************************
1093 : : // Decide if physics time falls into a field output time range
1094 : : //! \return True if physics time falls into a field output time range
1095 : : // *****************************************************************************
1096 : : {
1097 [ + + ]: 38466 : if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;
1098 : :
1099 : 32050 : const auto eps = std::numeric_limits< tk::real >::epsilon();
1100 : :
1101 : 32050 : bool output = false;
1102 : :
1103 : 32050 : const auto& rf = g_inputdeck.get< tag::field_output, tag::time_range >();
1104 [ - + ]: 32050 : if (!rf.empty()) {
1105 [ - - ][ - - ]: 0 : if (m_t > rf[0] and m_t < rf[1])
[ - - ]
1106 : 0 : output |= std::floor(m_t/rf[2]) - m_rangeFieldFloor > eps;
1107 : : }
1108 : :
1109 : 32050 : return output;
1110 : : }
1111 : :
1112 : : bool
1113 : 45134 : Discretization::histiter() const
1114 : : // *****************************************************************************
1115 : : // Decide if history output iteration count interval is hit
1116 : : //! \return True if history output iteration count interval is hit
1117 : : // *****************************************************************************
1118 : : {
1119 : 45134 : const auto hist = g_inputdeck.get< tag::history_output, tag::interval >();
1120 : 45134 : const auto& hist_points = g_inputdeck.get< tag::history_output, tag::point >();
1121 : :
1122 [ + + ][ + - ]: 45134 : return m_it % hist == 0 and not hist_points.empty();
1123 : : }
1124 : :
1125 : : bool
1126 : 44847 : Discretization::histtime() const
1127 : : // *****************************************************************************
1128 : : // Decide if history output physics time interval is hit
1129 : : //! \return True if history output physics time interval is hit
1130 : : // *****************************************************************************
1131 : : {
1132 [ + + ]: 44847 : if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;
1133 : :
1134 : 34927 : const auto eps = std::numeric_limits< tk::real >::epsilon();
1135 : 34927 : const auto ht = g_inputdeck.get< tag::history_output, tag::time_interval >();
1136 : :
1137 [ - + ]: 34927 : if (ht < eps) return false;
1138 : :
1139 : 34927 : return std::floor(m_t/ht) - m_physHistFloor > eps;
1140 : : }
1141 : :
1142 : : bool
1143 : 44841 : Discretization::histrange() const
1144 : : // *****************************************************************************
1145 : : // Decide if physics time falls into a history output time range
1146 : : //! \return True if physics time falls into a history output time range
1147 : : // *****************************************************************************
1148 : : {
1149 [ + + ]: 44841 : if (g_inputdeck.get< tag::cmd, tag::benchmark >()) return false;
1150 : :
1151 : 34921 : const auto eps = std::numeric_limits< tk::real >::epsilon();
1152 : :
1153 : 34921 : bool output = false;
1154 : :
1155 : 34921 : const auto& rh = g_inputdeck.get< tag::history_output, tag::time_range >();
1156 [ - + ]: 34921 : if (!rh.empty()) {
1157 [ - - ][ - - ]: 0 : if (m_t > rh[0] and m_t < rh[1])
[ - - ]
1158 : 0 : output |= std::floor(m_t/rh[2]) - m_rangeHistFloor > eps;
1159 : : }
1160 : :
1161 : 34921 : return output;
1162 : : }
1163 : :
1164 : : bool
1165 : 53928 : Discretization::finished() const
1166 : : // *****************************************************************************
1167 : : // Decide if this is the last time step
1168 : : //! \return True if this is the last time step
1169 : : // *****************************************************************************
1170 : : {
1171 : 53928 : const auto eps = std::numeric_limits< tk::real >::epsilon();
1172 : 53928 : const auto nstep = g_inputdeck.get< tag::nstep >();
1173 : 53928 : const auto term = g_inputdeck.get< tag::term >();
1174 : :
1175 [ + + ][ + + ]: 53928 : return std::abs(m_t-term) < eps or m_it >= nstep;
1176 : : }
1177 : :
1178 : : void
1179 : 42748 : Discretization::status()
1180 : : // *****************************************************************************
1181 : : // Output one-liner status report
1182 : : // *****************************************************************************
1183 : : {
1184 : : // Query after how many time steps user wants TTY dump
1185 : 42748 : const auto tty = g_inputdeck.get< tag::ttyi >();
1186 : :
1187 : : // estimate grind time (taken between this and the previous time step)
1188 : : using std::chrono::duration_cast;
1189 : : using ms = std::chrono::milliseconds;
1190 : : using clock = std::chrono::high_resolution_clock;
1191 [ + - ][ + - ]: 42748 : auto grind_time = duration_cast< ms >(clock::now() - m_prevstatus).count();
1192 : 42748 : m_prevstatus = clock::now();
1193 : :
1194 [ + + ][ + + ]: 42748 : if (thisIndex==0 and m_meshid == 0 and not (m_it%tty)) {
[ + + ]
1195 : :
1196 : 2386 : const auto term = g_inputdeck.get< tag::term >();
1197 : 2386 : const auto t0 = g_inputdeck.get< tag::t0 >();
1198 : 2386 : const auto nstep = g_inputdeck.get< tag::nstep >();
1199 : 2386 : const auto diag = g_inputdeck.get< tag::diagnostics, tag::interval >();
1200 : 2386 : const auto lbfreq = g_inputdeck.get< tag::cmd, tag::lbfreq >();
1201 : 2386 : const auto rsfreq = g_inputdeck.get< tag::cmd, tag::rsfreq >();
1202 : 2386 : const auto verbose = g_inputdeck.get< tag::cmd, tag::verbose >();
1203 : 2386 : const auto benchmark = g_inputdeck.get< tag::cmd, tag::benchmark >();
1204 : 2386 : const auto steady = g_inputdeck.get< tag::steady_state >();
1205 : :
1206 : : // estimate time elapsed and time for accomplishment
1207 [ + - ][ + - ]: 2386 : tk::Timer::Watch ete, eta;
1208 [ + + ][ + - ]: 2386 : if (not steady) m_timer.eta( term-t0, m_t-t0, nstep, m_it, ete, eta );
1209 : :
1210 : : const auto& def =
1211 : 2386 : g_inputdeck_defaults.get< tag::cmd, tag::io, tag::screen >();
1212 [ + - ]: 2386 : tk::Print print( g_inputdeck.get< tag::cmd >().logname( def, m_nrestart ),
1213 : : verbose ? std::cout : std::clog,
1214 [ + - ][ + - ]: 7158 : std::ios_base::app );
1215 : :
1216 : : // Output one-liner
1217 : 2386 : print << std::setfill(' ') << std::setw(8) << m_it << " "
1218 : 0 : << std::scientific << std::setprecision(6)
1219 : 0 : << std::setw(12) << m_t << " "
1220 : 2386 : << m_dt << " "
1221 : 0 : << std::setfill('0')
1222 : 0 : << std::setw(3) << ete.hrs.count() << ":"
1223 : 0 : << std::setw(2) << ete.min.count() << ":"
1224 : 0 : << std::setw(2) << ete.sec.count() << " "
1225 : 0 : << std::setw(3) << eta.hrs.count() << ":"
1226 : 0 : << std::setw(2) << eta.min.count() << ":"
1227 : 0 : << std::setw(2) << eta.sec.count() << " "
1228 : 0 : << std::scientific << std::setprecision(6) << std::setfill(' ')
1229 [ + - ][ + - ]: 2386 : << std::setw(9) << grind_time << " ";
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
1230 : :
1231 : : // Augment one-liner status with output indicators
1232 [ + - ][ + + ]: 2386 : if (fielditer() or fieldtime() or fieldrange()) print << 'f';
[ + - ][ + + ]
[ + - ][ - + ]
[ + + ][ + - ]
1233 [ + + ][ + - ]: 2386 : if (not (m_it % diag)) print << 'd';
1234 [ + - ][ + + ]: 2386 : if (histiter() or histtime() or histrange()) print << 't';
[ + - ][ + + ]
[ + - ][ - + ]
[ + + ][ + - ]
1235 [ - + ][ - - ]: 2386 : if (m_refined) print << 'h';
1236 [ + + ][ + - ]: 2386 : if (not (m_it % lbfreq) && not finished()) print << 'l';
[ + + ][ + + ]
[ + - ]
1237 [ + + ][ + - ]: 2386 : if (not benchmark && (not (m_it % rsfreq) || finished())) print << 'r';
[ + - ][ + + ]
[ + + ][ + - ]
1238 : :
1239 [ + + ][ + - ]: 2386 : if (not m_meshvel_converged) print << 'a';
1240 : 2386 : m_meshvel_converged = true; // get ready for next time step
1241 : :
1242 [ + - ]: 2386 : print << std::endl;
1243 : : }
1244 : 42748 : }
1245 : :
1246 : : #include "NoWarning/discretization.def.h"
|
