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