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 : : 2022-2024 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \details Data and functionality common to all discretization schemes
10 : : \see Discretization.h and Discretization.C for more info.
11 : : */
12 : : // *****************************************************************************
13 : :
14 : : #include <iomanip>
15 : :
16 : : #include "Reorder.hpp"
17 : : #include "Vector.hpp"
18 : : #include "DerivedData.hpp"
19 : : #include "Discretization.hpp"
20 : : #include "MeshWriter.hpp"
21 : : #include "DiagWriter.hpp"
22 : : #include "InciterConfig.hpp"
23 : : #include "Print.hpp"
24 : : #include "Around.hpp"
25 : : #include "XystBuildConfig.hpp"
26 : : #include "Box.hpp"
27 : :
28 : : namespace inciter {
29 : :
30 : : static CkReduction::reducerType PDFMerger;
31 : : extern ctr::Config g_cfg;
32 : :
33 : : } // inciter::
34 : :
35 : : using inciter::Discretization;
36 : :
37 : 2475 : Discretization::Discretization(
38 : : std::size_t meshid,
39 : : const CProxy_Transporter& transporter,
40 : : const tk::CProxy_MeshWriter& meshwriter,
41 : : const tk::UnsMesh::CoordMap& coordmap,
42 : : const tk::UnsMesh::Chunk& el,
43 : : const std::map< int, std::unordered_set< std::size_t > >& nodeCommMap,
44 : 2475 : int nc ) :
45 : 2475 : m_meshid( meshid ),
46 : 2475 : m_nchare( nc ),
47 : 2475 : m_it( 0 ),
48 : 2475 : m_itr( 0 ),
49 : 2475 : m_itf( 0 ),
50 : 2475 : m_initial( 1 ),
51 : 2475 : m_t( g_cfg.get< tag::t0 >() ),
52 : 2475 : m_lastDumpTime( -std::numeric_limits< tk::real >::max() ),
53 : 2475 : m_physFieldFloor( 0.0 ),
54 : 2475 : m_physHistFloor( 0.0 ),
55 : 2475 : m_physIntegFloor( 0.0 ),
56 [ + - ]: 2475 : m_rangeFieldFloor( g_cfg.get< tag::fieldout_range >().size(), 0.0 ),
57 [ + - ][ - - ]: 2475 : m_rangeHistFloor( g_cfg.get< tag::histout_range >().size(), 0.0 ),
58 [ + - ][ - - ]: 2475 : m_rangeIntegFloor( g_cfg.get< tag::integout_range >().size(), 0.0 ),
59 : 2475 : m_dt( g_cfg.get< tag::dt >() ),
60 : 2475 : m_dtn( m_dt ),
61 : 2475 : m_nvol( 0 ),
62 : : m_transporter( transporter ),
63 : : m_meshwriter( meshwriter ),
64 : : m_el( el ), // fills m_inpoel, m_gid, m_lid
65 : 2475 : m_coord( setCoord( coordmap ) ),
66 : 2475 : m_meshvol( 0.0 ),
67 [ + - ]: 2475 : m_v( m_gid.size(), 0.0 ),
68 [ + - ][ - - ]: 2475 : m_vol( m_gid.size(), 0.0 ),
69 : 2475 : m_refined( 0 ),
70 : 2475 : m_prevstatus( std::chrono::high_resolution_clock::now() ),
71 : 2475 : m_nrestart( 0 ),
72 : 2475 : m_res( 0.0 ),
73 : 2475 : m_res0( 0.0 ),
74 : 2475 : m_res1( 0.0 ),
75 : 2475 : m_dea( 0 ),
76 : 2475 : m_deastarted( 0 ),
77 : 2475 : m_pit( 0 ),
78 [ + - ][ + - ]: 12375 : m_mit( 0 )
[ + - ][ + - ]
79 : : // *****************************************************************************
80 : : // Constructor
81 : : //! \param[in] meshid Mesh ID
82 : : //! \param[in] transporter Host (Transporter) proxy
83 : : //! \param[in] meshwriter Mesh writer proxy
84 : : //! \param[in] coordmap Coordinates of mesh nodes and their global IDs
85 : : //! \param[in] el Elements of the mesh chunk we operate on
86 : : //! \param[in] nodeCommMap Node lists associated to chare IDs bordering the
87 : : //! mesh chunk we operate on
88 : : //! \param[in] nc Total number of Discretization chares
89 : : // *****************************************************************************
90 : : {
91 : : Assert( !m_inpoel.empty(), "No elements assigned to Discretization chare" );
92 : : Assert( tk::positiveJacobians( m_inpoel, m_coord ),
93 : : "Jacobian in input mesh to Discretization non-positive" );
94 : : #if not defined(__INTEL_COMPILER) || defined(NDEBUG)
95 : : // The above ifdef skips running the conformity test with the intel compiler
96 : : // in debug mode only. This is necessary because in tk::conforming(), filling
97 : : // up the map can fail with some meshes (only in parallel), e.g., tube.exo,
98 : : // used by some regression tests, due to the intel compiler generating some
99 : : // garbage incorrect code - only in debug, only in parallel, only with that
100 : : // mesh.
101 : : Assert( tk::conforming( m_inpoel, m_coord ),
102 : : "Input mesh to Discretization not conforming" );
103 : : #endif
104 : :
105 : : // Store node communication map
106 [ + + ]: 26921 : for (const auto& [c,map] : nodeCommMap) m_nodeCommMap[c] = map;
107 : :
108 : : // Get ready for computing/communicating nodal volumes
109 [ + - ]: 2475 : startvol();
110 : :
111 : : // Find host elements of user-specified points where time histories are
112 : : // saved, and save the shape functions evaluated at the point locations
113 : : const auto& pt = g_cfg.get< tag::histout >();
114 [ + + ]: 2715 : for (std::size_t p=0; p<pt.size(); ++p) {
115 : : std::array< tk::real, 4 > N;
116 : : const auto& l = pt[p];
117 [ + + ]: 117875 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
118 [ + - ][ + + ]: 117710 : if (tk::intet( m_coord, m_inpoel, l, e, N )) {
119 [ + - ]: 75 : m_histdata.push_back( HistData{{ "p"+std::to_string(p+1), e, N }} );
120 : 75 : break;
121 : : }
122 : : }
123 : : }
124 : :
125 : : // Compute number of mesh points owned
126 : 2475 : std::size_t npoin = m_gid.size();
127 [ + - ][ + + ]: 305305 : for (auto g : m_gid) if (tk::slave( m_nodeCommMap, g, thisIndex ) ) --npoin;
[ + + ]
128 : :
129 : : // Tell the RTS that the Discretization chares have been created and compute
130 : : // the total number of mesh points across the distributed mesh
131 [ + - ]: 2475 : std::vector< std::size_t > meshdata{ m_meshid, npoin };
132 [ + - ]: 2475 : contribute( meshdata, CkReduction::sum_ulong,
133 [ + - ][ - - ]: 2475 : CkCallback( CkReductionTarget(Transporter,disccreated), m_transporter ) );
134 : 2475 : }
135 : :
136 : : void
137 : 0 : Discretization::resizePostAMR(
138 : : const tk::UnsMesh::Chunk& chunk,
139 : : const tk::UnsMesh::Coords& coord,
140 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
141 : : nodeCommMap,
142 : : const std::set< std::size_t >& /*removedNodes*/ )
143 : : // *****************************************************************************
144 : : // Resize mesh data structures after mesh refinement
145 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
146 : : //! \param[in] coord New mesh node coordinates
147 : : //! \param[in] nodeCommMap New node communication map
148 : : //! \param[in] removedNodes Newly removed mesh node local ids
149 : : // *****************************************************************************
150 : : {
151 : : m_el = chunk; // updates m_inpoel, m_gid, m_lid
152 : : m_nodeCommMap = nodeCommMap;
153 : :
154 : : // Update mesh volume container size
155 : 0 : m_vol.resize( m_gid.size(), 0.0 );
156 : :
157 : : // update mesh node coordinates
158 : : m_coord = coord;
159 : 0 : }
160 : :
161 : : void
162 : 2475 : Discretization::startvol()
163 : : // *****************************************************************************
164 : : // Get ready for (re-)computing/communicating nodal volumes
165 : : // *****************************************************************************
166 : : {
167 : 2475 : m_nvol = 0;
168 [ + - ]: 4950 : thisProxy[ thisIndex ].wait4vol();
169 : :
170 : : // Zero out mesh volume container
171 : : std::fill( begin(m_vol), end(m_vol), 0.0 );
172 : :
173 : : // Clear receive buffer that will be used for collecting nodal volumes
174 : : m_volc.clear();
175 : 2475 : }
176 : :
177 : : void
178 : 783 : Discretization::registerReducers()
179 : : // *****************************************************************************
180 : : // Configure Charm++ reduction types
181 : : //! \details Since this is a [initnode] routine, see the .ci file, the
182 : : //! Charm++ runtime system executes the routine exactly once on every
183 : : //! logical node early on in the Charm++ init sequence. Must be static as
184 : : //! it is called without an object. See also: Section "Initializations at
185 : : //! Program Startup" at in the Charm++ manual
186 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
187 : : // *****************************************************************************
188 : : {
189 : 783 : PDFMerger = CkReduction::addReducer( tk::mergeUniPDFs );
190 : 783 : }
191 : :
192 : : tk::UnsMesh::Coords
193 : 2475 : Discretization::setCoord( const tk::UnsMesh::CoordMap& coordmap )
194 : : // *****************************************************************************
195 : : // Set mesh coordinates based on coordinates map
196 : : // *****************************************************************************
197 : : {
198 : : Assert( coordmap.size() == m_gid.size(), "Size mismatch" );
199 : : Assert( coordmap.size() == m_lid.size(), "Size mismatch" );
200 : :
201 : : tk::UnsMesh::Coords coord;
202 [ + - ]: 2475 : coord[0].resize( coordmap.size() );
203 [ + - ]: 2475 : coord[1].resize( coordmap.size() );
204 [ + - ]: 2475 : coord[2].resize( coordmap.size() );
205 : :
206 [ + + ]: 305305 : for (const auto& [ gid, coords ] : coordmap) {
207 : 302830 : auto i = tk::cref_find( m_lid, gid );
208 : 302830 : coord[0][i] = coords[0];
209 : 302830 : coord[1][i] = coords[1];
210 : 302830 : coord[2][i] = coords[2];
211 : : }
212 : :
213 : 2475 : return coord;
214 : 0 : }
215 : :
216 : : void
217 : 1581 : Discretization::remap(
218 : : const std::unordered_map< std::size_t, std::size_t >& map )
219 : : // *****************************************************************************
220 : : // Remap mesh data based on new local ids
221 : : //! \param[in] map Mapping of old->new local ids
222 : : // *****************************************************************************
223 : : {
224 : : // Remap connectivity containing local IDs
225 [ + + ]: 2646569 : for (auto& l : m_inpoel) l = tk::cref_find(map,l);
226 : :
227 : : // Remap global->local id map
228 [ + + ]: 187643 : for (auto& [g,l] : m_lid) l = tk::cref_find(map,l);
229 : :
230 : : // Remap global->local id map
231 : 1581 : auto maxid = std::numeric_limits< std::size_t >::max();
232 : 1581 : std::vector< std::size_t > newgid( m_gid.size(), maxid );
233 [ + + ]: 187643 : for (const auto& [o,n] : map) newgid[n] = m_gid[o];
234 : 1581 : m_gid = std::move( newgid );
235 : :
236 : : Assert( std::all_of( m_gid.cbegin(), m_gid.cend(),
237 : : [=](std::size_t i){ return i < maxid; } ),
238 : : "Not all gid have been remapped" );
239 : :
240 : : // Remap nodal volumes (with contributions along chare-boundaries)
241 [ + - ][ - - ]: 1581 : std::vector< tk::real > newvol( m_vol.size(), 0.0 );
242 [ + + ]: 187643 : for (const auto& [o,n] : map) newvol[n] = m_vol[o];
243 : 1581 : m_vol = std::move( newvol );
244 : :
245 : : // Remap nodal volumes (without contributions along chare-boundaries)
246 [ + - ][ - - ]: 1581 : std::vector< tk::real > newv( m_v.size(), 0.0 );
247 [ + + ]: 187643 : for (const auto& [o,n] : map) newv[n] = m_v[o];
248 : 1581 : m_v = std::move( newv );
249 : :
250 : : // Remap locations of node coordinates
251 : : tk::UnsMesh::Coords newcoord;
252 : : auto npoin = m_coord[0].size();
253 [ + - ]: 1581 : newcoord[0].resize( npoin );
254 [ + - ]: 1581 : newcoord[1].resize( npoin );
255 [ + - ]: 1581 : newcoord[2].resize( npoin );
256 [ + + ]: 187643 : for (const auto& [o,n] : map) {
257 : 186062 : newcoord[0][n] = m_coord[0][o];
258 : 186062 : newcoord[1][n] = m_coord[1][o];
259 : 186062 : newcoord[2][n] = m_coord[2][o];
260 : : }
261 : : m_coord = std::move( newcoord );
262 : 3162 : }
263 : :
264 : : void
265 : 2475 : Discretization::setRefiner( const CProxy_Refiner& ref )
266 : : // *****************************************************************************
267 : : // Set Refiner Charm++ proxy
268 : : //! \param[in] ref Incoming refiner proxy to store
269 : : // *****************************************************************************
270 : : {
271 : : m_refiner = ref;
272 : 2475 : }
273 : :
274 : : void
275 : 2475 : Discretization::vol()
276 : : // *****************************************************************************
277 : : // Sum mesh volumes to nodes, start communicating them on chare-boundaries
278 : : // *****************************************************************************
279 : : {
280 : : const auto& x = m_coord[0];
281 : : const auto& y = m_coord[1];
282 : : const auto& z = m_coord[2];
283 : :
284 : : // Compute nodal volumes on our chunk of the mesh
285 [ + + ]: 1102774 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
286 : 1100299 : const auto N = m_inpoel.data() + e*4;
287 : : const std::array< tk::real, 3 >
288 [ - + ]: 1100299 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
289 : 1100299 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
290 [ - + ]: 1100299 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
291 : 1100299 : const auto J = tk::triple( ba, ca, da ) / 24.0;
292 [ - + ][ - - ]: 1100299 : ErrChk( J > 0, "Element Jacobian non-positive: PE:" +
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
293 : : std::to_string(CkMyPe()) + ", node IDs: " +
294 : : std::to_string(m_gid[N[0]]) + ',' +
295 : : std::to_string(m_gid[N[1]]) + ',' +
296 : : std::to_string(m_gid[N[2]]) + ',' +
297 : : std::to_string(m_gid[N[3]]) + ", coords: (" +
298 : : std::to_string(x[N[0]]) + ", " +
299 : : std::to_string(y[N[0]]) + ", " +
300 : : std::to_string(z[N[0]]) + "), (" +
301 : : std::to_string(x[N[1]]) + ", " +
302 : : std::to_string(y[N[1]]) + ", " +
303 : : std::to_string(z[N[1]]) + "), (" +
304 : : std::to_string(x[N[2]]) + ", " +
305 : : std::to_string(y[N[2]]) + ", " +
306 : : std::to_string(z[N[2]]) + "), (" +
307 : : std::to_string(x[N[3]]) + ", " +
308 : : std::to_string(y[N[3]]) + ", " +
309 : : std::to_string(z[N[3]]) + ')' );
310 : : // scatter add V/4 to nodes
311 [ + + ]: 5501495 : for (std::size_t j=0; j<4; ++j) m_vol[N[j]] += J;
312 : : }
313 : :
314 : : // Store nodal volumes without contributions from other chares on
315 : : // chare-boundaries
316 : 2475 : m_v = m_vol;
317 : :
318 : : // Send our nodal volume contributions to neighbor chares
319 [ + + ]: 2475 : if (m_nodeCommMap.empty()) {
320 : 73 : comvol_complete();
321 : : } else {
322 [ + + ]: 26848 : for (const auto& [c,n] : m_nodeCommMap) {
323 : 24446 : std::vector< tk::real > v( n.size() );
324 : : std::size_t j = 0;
325 [ + + ]: 169856 : for (auto i : n) v[ j++ ] = m_vol[ tk::cref_find(m_lid,i) ];
326 [ + - ][ + - ]: 73338 : thisProxy[c].comvol( thisIndex,
[ + - ][ - - ]
327 [ + - ][ - + ]: 48892 : std::vector<std::size_t>(begin(n), end(n)), v );
[ - - ]
328 : : }
329 : : }
330 : :
331 : 2475 : ownvol_complete();
332 : 2475 : }
333 : :
334 : : void
335 : 24446 : Discretization::comvol( int c,
336 : : const std::vector< std::size_t >& gid,
337 : : const std::vector< tk::real >& nodevol )
338 : : // *****************************************************************************
339 : : // Receive nodal volumes on chare-boundaries
340 : : //! \param[in] c Sender chare id
341 : : //! \param[in] gid Global mesh node IDs at which we receive volume contributions
342 : : //! \param[in] nodevol Partial sums of nodal volume contributions to
343 : : //! chare-boundary nodes
344 : : // *****************************************************************************
345 : : {
346 : : Assert( nodevol.size() == gid.size(), "Size mismatch" );
347 : :
348 : : auto& cvolc = m_cvolc[c];
349 [ + + ]: 169856 : for (std::size_t i=0; i<gid.size(); ++i) {
350 : 145410 : m_volc[ gid[i] ] += nodevol[i];
351 : 145410 : cvolc[ tk::cref_find(m_lid,gid[i]) ] = nodevol[i];
352 : : }
353 : :
354 [ + + ]: 24446 : if (++m_nvol == m_nodeCommMap.size()) {
355 : 2402 : m_nvol = 0;
356 : 2402 : comvol_complete();
357 : : }
358 : 24446 : }
359 : :
360 : : void
361 : 2475 : Discretization::totalvol()
362 : : // *****************************************************************************
363 : : // Sum mesh volumes and contribute own mesh volume to total volume
364 : : // *****************************************************************************
365 : : {
366 : : // Add received contributions to nodal volumes
367 [ + + ]: 88193 : for (const auto& [gid, vol] : m_volc)
368 : 85718 : m_vol[ tk::cref_find(m_lid,gid) ] += vol;
369 : :
370 : : // Clear receive buffer
371 : 2475 : tk::destroy(m_volc);
372 : :
373 : : // Sum mesh volume to host
374 : : std::vector< tk::real > tvol{ 0.0,
375 : 2475 : static_cast<tk::real>(m_initial),
376 : 2475 : static_cast<tk::real>(m_meshid) };
377 [ + + ]: 305305 : for (auto v : m_v) tvol[0] += v;
378 [ + - ]: 2475 : contribute( tvol, CkReduction::sum_double,
379 [ + - ][ - - ]: 2475 : CkCallback(CkReductionTarget(Transporter,totalvol), m_transporter) );
380 : 2475 : }
381 : :
382 : : void
383 : 2475 : Discretization::stat( tk::real mesh_volume )
384 : : // *****************************************************************************
385 : : // Compute mesh cell statistics
386 : : //! \param[in] mesh_volume Total mesh volume
387 : : // *****************************************************************************
388 : : {
389 : : // Store total mesh volume
390 : 2475 : m_meshvol = mesh_volume;
391 : :
392 : : const auto& x = m_coord[0];
393 : : const auto& y = m_coord[1];
394 : : const auto& z = m_coord[2];
395 : :
396 : 2475 : auto MIN = -std::numeric_limits< tk::real >::max();
397 : 2475 : auto MAX = std::numeric_limits< tk::real >::max();
398 : 2475 : std::vector< tk::real > min( 6, MAX );
399 [ + - ][ - - ]: 2475 : std::vector< tk::real > max( 6, MIN );
400 [ + - ][ + - ]: 2475 : std::vector< tk::real > sum( 9, 0.0 );
[ - - ]
401 : : tk::UniPDF edgePDF( 1e-4 );
402 : : tk::UniPDF volPDF( 1e-4 );
403 : : tk::UniPDF ntetPDF( 1e-4 );
404 : :
405 : : // Compute points surrounding points
406 [ + - ][ + - ]: 2475 : auto psup = tk::genPsup( m_inpoel, 4, tk::genEsup(m_inpoel,4) );
407 : : Assert( psup.second.size()-1 == m_gid.size(),
408 : : "Number of mesh points and number of global IDs unequal" );
409 : :
410 : : // Compute edge length statistics
411 : : // Note that while the min and max edge lengths are independent of the number
412 : : // of CPUs (by the time they are aggregated across all chares), the sum of
413 : : // the edge lengths and the edge length PDF are not. This is because the
414 : : // edges on the chare-boundary are counted multiple times and we
415 : : // conscientiously do not make an effort to precisely compute this, because
416 : : // that would require communication and more complex logic. Since these
417 : : // statistics are intended as simple average diagnostics, we ignore these
418 : : // small differences. For reproducible average edge lengths and edge length
419 : : // PDFs, run the mesh in serial.
420 : : tk::UnsMesh::EdgeSet edges;
421 [ + + ]: 305305 : for (std::size_t p=0; p<m_gid.size(); ++p) {
422 [ + + ][ + + ]: 3485478 : for (auto i : tk::Around(psup,p)) {
423 : 3182648 : const auto dx = x[ i ] - x[ p ];
424 : 3182648 : const auto dy = y[ i ] - y[ p ];
425 : 3182648 : const auto dz = z[ i ] - z[ p ];
426 [ + + ]: 3182648 : const auto length = std::sqrt( dx*dx + dy*dy + dz*dz );
427 [ + + ]: 3182648 : if (length < min[0]) min[0] = length;
428 [ + + ]: 3182648 : if (length > max[0]) max[0] = length;
429 : 3182648 : sum[0] += 1.0;
430 : 3182648 : sum[1] += length;
431 [ + - ]: 3182648 : edgePDF.add( length );
432 [ + - ]: 3182648 : edges.insert( { m_gid[i], m_gid[p] } );
433 : : }
434 : : }
435 : :
436 : : // Compute mesh cell volume statistics
437 [ + + ]: 1102774 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
438 [ + + ]: 1100299 : const std::array< std::size_t, 4 > N{{ m_inpoel[e*4+0], m_inpoel[e*4+1],
439 [ + + ]: 1100299 : m_inpoel[e*4+2], m_inpoel[e*4+3] }};
440 : : const std::array< tk::real, 3 >
441 : 1100299 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
442 : 1100299 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
443 [ + + ]: 1100299 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
444 : 1100299 : const auto L = std::cbrt( tk::triple( ba, ca, da ) / 6.0 );
445 [ + + ]: 1100299 : if (L < min[1]) min[1] = L;
446 [ + + ]: 1100299 : if (L > max[1]) max[1] = L;
447 : 1100299 : sum[2] += 1.0;
448 : 1100299 : sum[3] += L;
449 [ + - ]: 1100299 : volPDF.add( L );
450 : : }
451 : :
452 : : // Contribute statistics
453 : 2475 : sum[4] = 1.0;
454 : 2475 : min[2] = max[2] = sum[5] = static_cast< tk::real >( m_inpoel.size() / 4 );
455 : 2475 : min[3] = max[3] = sum[6] = static_cast< tk::real >( m_gid.size() );
456 : 2475 : min[4] = max[4] = sum[7] = static_cast< tk::real >( edges.size() );
457 : 2475 : min[5] = max[5] = sum[8] =
458 : 2475 : static_cast< tk::real >( tk::sumvalsize(m_nodeCommMap) ) /
459 : : static_cast< tk::real >( m_gid.size() );
460 [ + - ]: 2475 : ntetPDF.add( min[2] );
461 : :
462 [ + - ]: 2475 : min.push_back( static_cast<tk::real>(m_meshid) );
463 [ + - ]: 2475 : max.push_back( static_cast<tk::real>(m_meshid) );
464 [ + - ]: 2475 : sum.push_back( static_cast<tk::real>(m_meshid) );
465 : :
466 : : // Contribute to mesh statistics across all Discretization chares
467 [ + - ]: 2475 : contribute( min, CkReduction::min_double,
468 [ + - ]: 2475 : CkCallback(CkReductionTarget(Transporter,minstat), m_transporter) );
469 : : contribute( max, CkReduction::max_double,
470 [ + - ]: 2475 : CkCallback(CkReductionTarget(Transporter,maxstat), m_transporter) );
471 : : contribute( sum, CkReduction::sum_double,
472 : 2475 : CkCallback(CkReductionTarget(Transporter,sumstat), m_transporter) );
473 : :
474 : : // Serialize PDFs to raw stream
475 [ + - ][ + - ]: 9900 : auto stream = tk::serialize( m_meshid, { edgePDF, volPDF, ntetPDF } );
[ + + ][ + - ]
[ - - ]
476 : : // Create Charm++ callback function for reduction of PDFs with
477 : : // Transporter::pdfstat() as the final target where the results will appear.
478 : : CkCallback cb( CkIndex_Transporter::pdfstat(nullptr), m_transporter );
479 : : // Contribute serialized PDF of partial sums to host via Charm++ reduction
480 [ + - ]: 2475 : contribute( stream.first, stream.second.get(), PDFMerger, cb );
481 [ + - ][ + - ]: 7425 : }
[ + - ][ - - ]
482 : :
483 : : void
484 : 2475 : Discretization::boxvol()
485 : : // *****************************************************************************
486 : : // Compute total box IC volume
487 : : // *****************************************************************************
488 : : {
489 : : // Determine which nodes reside in user-defined IC box(es) if any
490 : 2475 : m_boxnodes = problems::boxnodes( m_coord );
491 : :
492 : : // Compute partial box IC volume (just add up all boxes)
493 : : tk::real boxvol = 0.0;
494 : : // cppcheck-suppress useStlAlgorithm
495 [ + + ][ + + ]: 3954 : for (const auto& b : m_boxnodes) for (auto i : b) boxvol += m_v[i];
496 : :
497 : : // Sum up box IC volume across all chares
498 : 2475 : std::vector< tk::real > meshdata{ boxvol, static_cast<tk::real>(m_meshid) };
499 [ + - ]: 2475 : contribute( meshdata, CkReduction::sum_double,
500 [ + - ][ - - ]: 2475 : CkCallback(CkReductionTarget(Transporter,boxvol), m_transporter) );
501 : 2475 : }
502 : :
503 : : void
504 : 3347 : Discretization::write(
505 : : const std::vector< std::size_t >& inpoel,
506 : : const tk::UnsMesh::Coords& coord,
507 : : const std::map< int, std::vector< std::size_t > >& bface,
508 : : const std::map< int, std::vector< std::size_t > >& bnode,
509 : : const std::vector< std::size_t >& triinpoel,
510 : : const std::vector< std::string>& elemfieldnames,
511 : : const std::vector< std::string>& nodefieldnames,
512 : : const std::vector< std::string>& elemsurfnames,
513 : : const std::vector< std::string>& nodesurfnames,
514 : : const std::vector< std::vector< tk::real > >& elemfields,
515 : : const std::vector< std::vector< tk::real > >& nodefields,
516 : : const std::vector< std::vector< tk::real > >& elemsurfs,
517 : : const std::vector< std::vector< tk::real > >& nodesurfs,
518 : : CkCallback c )
519 : : // *****************************************************************************
520 : : // Output mesh and fields data (solution dump) to file(s)
521 : : //! \param[in] inpoel Mesh connectivity for the mesh chunk to be written
522 : : //! \param[in] coord Node coordinates of the mesh chunk to be written
523 : : //! \param[in] bface Map of boundary-face lists mapped to corresponding side set
524 : : //! ids for this mesh chunk
525 : : //! \param[in] bnode Map of boundary-node lists mapped to corresponding side set
526 : : //! ids for this mesh chunk
527 : : //! \param[in] triinpoel Interconnectivity of points and boundary-face in this
528 : : //! mesh chunk
529 : : //! \param[in] elemfieldnames Names of element fields to be output to file
530 : : //! \param[in] nodefieldnames Names of node fields to be output to file
531 : : //! \param[in] elemsurfnames Names of elem surface fields to be output to file
532 : : //! \param[in] nodesurfnames Names of node surface fields to be output to file
533 : : //! \param[in] elemfields Field data in mesh elements to output to file
534 : : //! \param[in] nodefields Field data in mesh nodes to output to file
535 : : //! \param[in] elemsurfs Surface field data in mesh elements to output to file
536 : : //! \param[in] nodesurfs Surface field data in mesh nodes to output to file
537 : : //! \param[in] c Function to continue with after the write
538 : : //! \details Since m_meshwriter is a Charm++ chare group, it never migrates and
539 : : //! an instance is guaranteed on every PE. We index the first PE on every
540 : : //! logical compute node. In Charm++'s non-SMP mode, a node is the same as a
541 : : //! PE, so the index is the same as CkMyPe(). In SMP mode the index is the
542 : : //! first PE on every logical node. In non-SMP mode this yields one or more
543 : : //! output files per PE with zero or non-zero virtualization, respectively. If
544 : : //! there are multiple chares on a PE, the writes are serialized per PE, since
545 : : //! only a single entry method call can be executed at any given time. In SMP
546 : : //! mode, still the same number of files are output (one per chare), but the
547 : : //! output is serialized through the first PE of each compute node. In SMP
548 : : //! mode, channeling multiple files via a single PE on each node is required
549 : : //! by NetCDF and HDF5, as well as ExodusII, since none of these libraries are
550 : : //! thread-safe.
551 : : // *****************************************************************************
552 : : {
553 : : // If the previous iteration refined (or moved) the mesh or this is called
554 : : // before the first time step, we also output the mesh.
555 : 3347 : bool meshoutput = m_itf == 0 ? true : false;
556 : :
557 : : auto eps = std::numeric_limits< tk::real >::epsilon();
558 : 3347 : bool fieldoutput = false;
559 : :
560 : : // Output field data only if there is no dump at this physical time yet
561 [ + - ]: 3347 : if (std::abs(m_lastDumpTime - m_t) > eps ) {
562 : 3347 : m_lastDumpTime = m_t;
563 : 3347 : ++m_itf;
564 : 3347 : fieldoutput = true;
565 : : }
566 : :
567 : : const auto& f = g_cfg.get< tag::fieldout >();
568 : 3347 : std::set< int > outsets( begin(f), end(f) );
569 : :
570 : : m_meshwriter[ CkNodeFirst( CkMyNode() ) ].
571 [ + - ][ + - ]: 10041 : write( m_meshid, meshoutput, fieldoutput, m_itr, m_itf, m_t, thisIndex,
572 : : g_cfg.get< tag::output >(),
573 : : inpoel, coord, bface, bnode, triinpoel, elemfieldnames,
574 : : nodefieldnames, elemsurfnames, nodesurfnames, elemfields, nodefields,
575 : : elemsurfs, nodesurfs, outsets, c );
576 : 3347 : }
577 : :
578 : : void
579 : 38556 : Discretization::setdt( tk::real newdt )
580 : : // *****************************************************************************
581 : : // Set time step size
582 : : //! \param[in] newdt Size of the new time step
583 : : // *****************************************************************************
584 : : {
585 : 38556 : m_dtn = m_dt;
586 : 38556 : m_dt = newdt;
587 : :
588 : : // Truncate the size of last time step
589 : 38556 : const auto term = g_cfg.get< tag::term >();
590 [ + + ]: 38556 : if (m_t+m_dt > term) m_dt = term - m_t;
591 : 38556 : }
592 : :
593 : : void
594 : 38588 : Discretization::next()
595 : : // *****************************************************************************
596 : : // Prepare for next step
597 : : // *****************************************************************************
598 : : {
599 : : // Update floor of physics time divided by output interval times
600 : : const auto eps = std::numeric_limits< tk::real >::epsilon();
601 : 38588 : const auto ft = g_cfg.get< tag::fieldout_time >();
602 [ + - ]: 38588 : if (ft > eps) m_physFieldFloor = std::floor( m_t / ft );
603 : 38588 : const auto ht = g_cfg.get< tag::histout_time >();
604 [ + - ]: 38588 : if (ht > eps) m_physHistFloor = std::floor( m_t / ht );
605 : 38588 : const auto it = g_cfg.get< tag::integout_time >();
606 [ + - ]: 38588 : if (it > eps) m_physIntegFloor = std::floor( m_t / it );
607 : :
608 : : // Update floors of physics time divided by output interval times for ranges
609 : : const auto& rf = g_cfg.get< tag::fieldout_range >();
610 [ + + ]: 38988 : for (std::size_t i=0; i<rf.size(); ++i) {
611 [ + + ][ + + ]: 400 : if (m_t > rf[i][0] and m_t < rf[i][1]) {
612 : 32 : m_rangeFieldFloor[i] = std::floor( m_t / rf[i][2] );
613 : : }
614 : : }
615 : :
616 : : const auto& rh = g_cfg.get< tag::histout_range >();
617 [ + + ]: 40828 : for (std::size_t i=0; i<rh.size(); ++i) {
618 [ + + ][ + + ]: 2240 : if (m_t > rh[i][0] and m_t < rh[i][1]) {
619 : 620 : m_rangeHistFloor[i] = std::floor( m_t / rh[i][2] );
620 : : }
621 : : }
622 : :
623 : : const auto& ri = g_cfg.get< tag::integout_range >();
624 [ + + ]: 38628 : for (std::size_t i=0; i<ri.size(); ++i) {
625 [ + + ][ + + ]: 40 : if (m_t > ri[i][0] and m_t < ri[i][1]) {
626 : 10 : m_rangeIntegFloor[i] = std::floor( m_t / ri[i][2] );
627 : : }
628 : : }
629 : :
630 : 38588 : ++m_it;
631 : 38588 : m_t += m_dt;
632 : 38588 : }
633 : :
634 : : void
635 : 2473 : Discretization::grindZero()
636 : : // *****************************************************************************
637 : : // Zero grind-time
638 : : // *****************************************************************************
639 : : {
640 : 2473 : m_prevstatus = std::chrono::high_resolution_clock::now();
641 : :
642 [ + + ][ + - ]: 2473 : if (thisIndex == 0 && m_meshid == 0) {
643 : : tk::Print() << "Starting time stepping ...\n";
644 : : }
645 : 2473 : }
646 : :
647 : : bool
648 : 36113 : Discretization::restarted( int nrestart )
649 : : // *****************************************************************************
650 : : // Detect if just returned from a checkpoint and if so, zero timers
651 : : //! \param[in] nrestart Number of times restarted
652 : : //! \return True if restart detected
653 : : // *****************************************************************************
654 : : {
655 : : // Detect if just restarted from checkpoint:
656 : : // nrestart == -1 if there was no checkpoint this step
657 : : // m_nrestart == nrestart if there was a checkpoint this step
658 : : // if both false, just restarted from a checkpoint
659 [ + + ][ + - ]: 36113 : bool restarted = nrestart != -1 and m_nrestart != nrestart;
660 : :
661 : : // If just restarted from checkpoint
662 : : if (restarted) {
663 : : // Update number of restarts
664 : 30 : m_nrestart = nrestart;
665 : : // Start timer measuring time stepping wall clock time
666 : 30 : m_timer.zero();
667 : : // Zero grind-timer
668 : 30 : grindZero();
669 : : }
670 : :
671 : 36113 : return restarted;
672 : : }
673 : :
674 : : std::string
675 : 726 : Discretization::histfilename( const std::string& id, std::streamsize precision )
676 : : // *****************************************************************************
677 : : // Construct history output filename
678 : : //! \param[in] id History point id
679 : : //! \param[in] precision Floating point precision to use for output
680 : : //! \return History file name
681 : : // *****************************************************************************
682 : : {
683 : : auto of = g_cfg.get< tag::output >();
684 [ + - ]: 726 : std::stringstream ss;
685 : :
686 [ + - ]: 726 : ss << std::setprecision(static_cast<int>(precision)) << of << ".hist." << id;
687 : :
688 : 726 : return ss.str();
689 : 726 : }
690 : :
691 : : void
692 : 74 : Discretization::histheader( std::vector< std::string >&& names )
693 : : // *****************************************************************************
694 : : // Output headers for time history files (one for each point)
695 : : //! \param[in] names History output variable names
696 : : // *****************************************************************************
697 : : {
698 [ + + ]: 149 : for (const auto& h : m_histdata) {
699 : 75 : auto prec = g_cfg.get< tag::histout_precision >();
700 : 75 : tk::DiagWriter hw( histfilename( h.get< tag::id >(), prec ),
701 : : g_cfg.get< tag::histout_format >(),
702 [ + - ]: 75 : prec );
703 [ + - ]: 75 : hw.header( names );
704 : : }
705 : 74 : }
706 : :
707 : : void
708 : 1238 : Discretization::history( std::vector< std::vector< tk::real > >&& data )
709 : : // *****************************************************************************
710 : : // Output time history for a time step
711 : : //! \param[in] data Time history data for all variables and equations integrated
712 : : // *****************************************************************************
713 : : {
714 : : Assert( data.size() == m_histdata.size(), "Size mismatch" );
715 : :
716 : : std::size_t i = 0;
717 [ + + ]: 1889 : for (const auto& h : m_histdata) {
718 : 651 : auto prec = g_cfg.get< tag::histout_precision >();
719 [ + - ]: 1302 : tk::DiagWriter hw( histfilename( h.get< tag::id >(), prec ),
720 : : g_cfg.get< tag::histout_format >(),
721 : : prec,
722 [ + - ]: 651 : std::ios_base::app );
723 [ + - ]: 651 : hw.write( m_it, m_t, m_dt, data[i] );
724 : 651 : ++i;
725 : : }
726 : 1238 : }
727 : :
728 : : bool
729 : 41468 : Discretization::fielditer() const
730 : : // *****************************************************************************
731 : : // Decide if field output iteration count interval is hit
732 : : //! \return True if field output iteration count interval is hit
733 : : // *****************************************************************************
734 : : {
735 [ + + ]: 41468 : if (g_cfg.get< tag::benchmark >()) return false;
736 : :
737 : 22868 : return m_it % g_cfg.get< tag::fieldout_iter >() == 0;
738 : : }
739 : :
740 : : bool
741 : 38877 : Discretization::fieldtime() const
742 : : // *****************************************************************************
743 : : // Decide if field output physics time interval is hit
744 : : //! \return True if field output physics time interval is hit
745 : : // *****************************************************************************
746 : : {
747 [ + + ]: 38877 : if (g_cfg.get< tag::benchmark >()) return false;
748 : :
749 : : auto eps = std::numeric_limits< tk::real >::epsilon();
750 : 20277 : auto ft = g_cfg.get< tag::fieldout_time >();
751 : :
752 [ + - ]: 20277 : if (ft < eps) return false;
753 : :
754 : 20277 : return std::floor(m_t/ft) - m_physFieldFloor > eps;
755 : : }
756 : :
757 : : bool
758 : 38853 : Discretization::fieldrange() const
759 : : // *****************************************************************************
760 : : // Decide if physics time falls into a field output time range
761 : : //! \return True if physics time falls into a field output time range
762 : : // *****************************************************************************
763 : : {
764 [ + + ]: 38853 : if (g_cfg.get< tag::benchmark >()) return false;
765 : :
766 : : const auto eps = std::numeric_limits< tk::real >::epsilon();
767 : :
768 : : bool output = false;
769 : :
770 : : const auto& rf = g_cfg.get< tag::fieldout_range >();
771 [ + + ]: 20753 : for (std::size_t i=0; i<rf.size(); ++i) {
772 [ + + ][ + + ]: 500 : if (m_t > rf[i][0] and m_t < rf[i][1]) {
773 : 40 : output |= std::floor(m_t/rf[i][2]) - m_rangeFieldFloor[i] > eps;
774 : : }
775 : : }
776 : :
777 : : return output;
778 : : }
779 : :
780 : : bool
781 : 41468 : Discretization::histiter() const
782 : : // *****************************************************************************
783 : : // Decide if history output iteration count interval is hit
784 : : //! \return True if history output iteration count interval is hit
785 : : // *****************************************************************************
786 : : {
787 [ + + ]: 41468 : if (g_cfg.get< tag::benchmark >()) return false;
788 : :
789 : 22868 : auto hist = g_cfg.get< tag::histout_iter >();
790 : : const auto& hist_points = g_cfg.get< tag::histout >();
791 : :
792 [ + + ][ - + ]: 22868 : return m_it % hist == 0 and not hist_points.empty();
793 : : }
794 : :
795 : : bool
796 : 40682 : Discretization::histtime() const
797 : : // *****************************************************************************
798 : : // Decide if history output physics time interval is hit
799 : : //! \return True if history output physics time interval is hit
800 : : // *****************************************************************************
801 : : {
802 [ + + ]: 40682 : if (g_cfg.get< tag::benchmark >()) return false;
803 : :
804 : : auto eps = std::numeric_limits< tk::real >::epsilon();
805 : 22082 : auto ht = g_cfg.get< tag::histout_time >();
806 : :
807 [ + - ]: 22082 : if (ht < eps) return false;
808 : :
809 : 22082 : return std::floor(m_t/ht) - m_physHistFloor > eps;
810 : : }
811 : :
812 : : bool
813 : 40424 : Discretization::histrange() const
814 : : // *****************************************************************************
815 : : // Decide if physics time falls into a history output time range
816 : : //! \return True if physics time falls into a history output time range
817 : : // *****************************************************************************
818 : : {
819 [ + + ]: 40424 : if (g_cfg.get< tag::benchmark >()) return false;
820 : :
821 : : auto eps = std::numeric_limits< tk::real >::epsilon();
822 : :
823 : : bool output = false;
824 : :
825 : : const auto& rh = g_cfg.get< tag::histout_range >();
826 [ + + ]: 23984 : for (std::size_t i=0; i<rh.size(); ++i) {
827 [ + + ][ + + ]: 2160 : if (m_t > rh[i][0] and m_t < rh[i][1]) {
828 : 590 : output |= std::floor(m_t/rh[i][2]) - m_rangeHistFloor[i] > eps;
829 : : }
830 : : }
831 : :
832 : : return output;
833 : : }
834 : :
835 : : bool
836 : 41468 : Discretization::integiter() const
837 : : // *****************************************************************************
838 : : // Decide if integral output iteration count interval is hit
839 : : //! \return True if integral output iteration count interval is hit
840 : : // *****************************************************************************
841 : : {
842 [ + + ]: 41468 : if (g_cfg.get< tag::benchmark >()) return false;
843 : :
844 : 22868 : auto integ = g_cfg.get< tag::integout_iter >();
845 : : const auto& sidesets_integral = g_cfg.get< tag::integout >();
846 : :
847 [ + + ][ - + ]: 22868 : return m_it % integ == 0 and not sidesets_integral.empty();
848 : : }
849 : :
850 : : bool
851 : 41182 : Discretization::integtime() const
852 : : // *****************************************************************************
853 : : // Decide if integral output physics time interval is hit
854 : : //! \return True if integral output physics time interval is hit
855 : : // *****************************************************************************
856 : : {
857 [ + + ]: 41182 : if (g_cfg.get< tag::benchmark >()) return false;
858 : :
859 : : auto eps = std::numeric_limits< tk::real >::epsilon();
860 : 22582 : auto it = g_cfg.get< tag::integout_time >();
861 : :
862 [ + - ]: 22582 : if (it < eps) return false;
863 : :
864 : 22582 : return std::floor(m_t/it) - m_physIntegFloor > eps;
865 : : }
866 : :
867 : : bool
868 : 41155 : Discretization::integrange() const
869 : : // *****************************************************************************
870 : : // Decide if physics time falls into a integral output time range
871 : : //! \return True if physics time falls into a integral output time range
872 : : // *****************************************************************************
873 : : {
874 [ + + ]: 41155 : if (g_cfg.get< tag::benchmark >()) return false;
875 : :
876 : : auto eps = std::numeric_limits< tk::real >::epsilon();
877 : :
878 : : bool output = false;
879 : :
880 : : const auto& ri = g_cfg.get< tag::integout_range >();
881 [ + + ]: 22615 : for (std::size_t i=0; i<ri.size(); ++i) {
882 [ + + ][ + + ]: 60 : if (m_t > ri[i][0] and m_t < ri[i][1]) {
883 : 15 : output |= std::floor(m_t/ri[i][2]) - m_rangeIntegFloor[i] > eps;
884 : : }
885 : : }
886 : :
887 : : return output;
888 : : }
889 : :
890 : : bool
891 : 43143 : Discretization::finished() const
892 : : // *****************************************************************************
893 : : // Decide if this is the last time step
894 : : //! \return True if this is the last time step
895 : : // *****************************************************************************
896 : : {
897 : : auto eps = std::numeric_limits< tk::real >::epsilon();
898 : 43143 : auto nstep = g_cfg.get< tag::nstep >();
899 : 43143 : auto term = g_cfg.get< tag::term >();
900 : 43143 : auto residual = g_cfg.get< tag::residual >();
901 : :
902 [ + + ][ + + ]: 43143 : return std::abs(m_t-term) < eps or m_it >= nstep or
903 [ + + ][ - + ]: 40263 : (m_res > 0.0 and m_res < residual);
904 : : }
905 : :
906 : : void
907 : 1060 : Discretization::residual( tk::real r )
908 : : // *****************************************************************************
909 : : // Update residual (during convergence to steady state)
910 : : //! \param[in] r Current residual
911 : : // *****************************************************************************
912 : : {
913 : 1060 : auto ttyi = g_cfg.get< tag::ttyi >();
914 : :
915 [ + - ]: 1060 : if (m_it % ttyi == 0) {
916 : 1060 : m_res0 = m_res1;
917 : 1060 : m_res1 = r;
918 : : }
919 : :
920 : 1060 : m_res = r;
921 : 1060 : }
922 : :
923 : : bool
924 : 55 : Discretization::deastart()
925 : : // *****************************************************************************
926 : : // Decide whether to start the deactivation procedure in this time step
927 : : //! \return True to start the deactivation prcedure in this time step
928 : : // *****************************************************************************
929 : : {
930 [ + + ][ + - ]: 55 : if (not m_deastarted and m_t > g_cfg.get<tag::deatime>() + m_dt) {
931 : 19 : m_deastarted = 1;
932 : 19 : return true;
933 : : }
934 : :
935 : : return false;
936 : : }
937 : :
938 : : bool
939 : 6086 : Discretization::deactivate() const
940 : : // *****************************************************************************
941 : : // Decide whether to run deactivation procedure in this time step
942 : : //! \return True to run deactivation prcedure in this time step
943 : : // *****************************************************************************
944 : : {
945 : 6086 : auto dea = g_cfg.get< tag::deactivate >();
946 : 6086 : auto deafreq = g_cfg.get< tag::deafreq >();
947 : :
948 [ + + ][ + - ]: 6086 : if (dea and !m_nodeCommMap.empty() and m_it % deafreq == 0)
[ + + ]
949 : : return true;
950 : : else
951 : 5886 : return false;
952 : : }
953 : :
954 : : void
955 : 10 : Discretization::deactivated( int n )
956 : : // *****************************************************************************
957 : : // Receive deactivation report
958 : : //! \param[in] n Sum of deactivated chares
959 : : // *****************************************************************************
960 : : {
961 : 10 : m_dea = n;
962 : 10 : }
963 : :
964 : : bool
965 : 17032 : Discretization::lb() const
966 : : // *****************************************************************************
967 : : // Decide whether to do load balancing this time step
968 : : //! \return True to do load balancing in this time step
969 : : // *****************************************************************************
970 : : {
971 : 17032 : auto lbfreq = g_cfg.get< tag::lbfreq >();
972 : 17032 : auto lbtime = g_cfg.get< tag::lbtime >();
973 : :
974 [ + - ][ + + ]: 17032 : if ((m_t > lbtime and m_it % lbfreq == 0) or m_it == 2)
[ + + ]
975 : : return true;
976 : : else
977 : 5114 : return false;
978 : : }
979 : :
980 : : void
981 : 563 : Discretization::pit( std::size_t it )
982 : : // *****************************************************************************
983 : : // Update number of pressure linear solve iterations taken
984 : : //! \param[in] it Number of pressure linear solve iterations taken
985 : : // *****************************************************************************
986 : : {
987 : 563 : m_pit = it;
988 : 563 : }
989 : :
990 : : void
991 : 20 : Discretization::mit( std::size_t it )
992 : : // *****************************************************************************
993 : : // Update number of momentum/transport linear solve iterations taken
994 : : //! \param[in] it Number of momentum/transport linear solve iterations taken
995 : : // *****************************************************************************
996 : : {
997 : 20 : m_mit = it;
998 : 20 : }
999 : :
1000 : : void
1001 : 246 : Discretization::npoin( std::size_t n )
1002 : : // *****************************************************************************
1003 : : // Set number of mesh points (across all meshes)
1004 : : //! \param[in] n Number of mesh points
1005 : : // *****************************************************************************
1006 : : {
1007 : 246 : m_npoin = n;
1008 : 246 : }
1009 : :
1010 : : void
1011 : 18397 : Discretization::status()
1012 : : // *****************************************************************************
1013 : : // Output one-liner status report
1014 : : // *****************************************************************************
1015 : : {
1016 : 18397 : const auto ttyi = g_cfg.get< tag::ttyi >();
1017 : :
1018 [ + + ][ + - ]: 18397 : if (thisIndex == 0 and m_meshid == 0 and m_it % ttyi == 0) {
[ + + ]
1019 : :
1020 : : // estimate grind time (taken between this and the previous time step)
1021 : : using std::chrono::duration_cast;
1022 : : using us = std::chrono::microseconds;
1023 : : using clock = std::chrono::high_resolution_clock;
1024 : 2880 : auto grind_time = duration_cast< us >(clock::now() - m_prevstatus).count()
1025 : 2880 : / static_cast< long >( ttyi );
1026 : 2880 : auto grind_perf = static_cast<tk::real>(grind_time)
1027 : 2880 : / static_cast<tk::real>(m_npoin);
1028 : 2880 : m_prevstatus = clock::now();
1029 : :
1030 : 2880 : const auto term = g_cfg.get< tag::term >();
1031 : 2880 : const auto t0 = g_cfg.get< tag::t0 >();
1032 : 2880 : const auto nstep = g_cfg.get< tag::nstep >();
1033 : 2880 : const auto diag = g_cfg.get< tag::diag_iter >();
1034 : 2880 : const auto rsfreq = g_cfg.get< tag::rsfreq >();
1035 : 2880 : const auto benchmark = g_cfg.get< tag::benchmark >();
1036 : 2880 : const auto residual = g_cfg.get< tag::residual >();
1037 : : const auto solver = g_cfg.get< tag::solver >();
1038 [ + + ]: 2880 : const auto pre = solver == "chocg" ? 1 : 0;
1039 : 2880 : const auto theta = g_cfg.get< tag::theta >();
1040 : : const auto eps = std::numeric_limits< tk::real >::epsilon();
1041 [ + + ][ + + ]: 2880 : const auto mom = solver == "chocg" and theta > eps ? 1 : 0;
1042 : :
1043 : : // estimate time elapsed and time for accomplishment
1044 : : tk::Timer::Watch ete, eta;
1045 [ + - ]: 2880 : m_timer.eta( term-t0, m_t-t0, nstep, m_it, m_res0, m_res1, residual,
1046 : : ete, eta );
1047 : :
1048 : : tk::Print print;
1049 : :
1050 : : // Output one-liner
1051 : 5760 : print << std::setfill(' ') << std::setw(8) << m_it << " "
1052 : : << std::scientific << std::setprecision(6)
1053 : : << std::setw(12) << m_t << " "
1054 : : << m_dt << " "
1055 : 5760 : << std::setfill('0')
1056 : : << std::setw(3) << ete.hrs.count() << ":"
1057 : : << std::setw(2) << ete.min.count() << ":"
1058 : : << std::setw(2) << ete.sec.count() << " "
1059 : : << std::setw(3) << eta.hrs.count() << ":"
1060 : : << std::setw(2) << eta.min.count() << ":"
1061 : : << std::setw(2) << eta.sec.count() << " "
1062 [ - - ]: 0 : << std::scientific << std::setprecision(6) << std::setfill(' ')
1063 : : << std::setw(9) << grind_time << " "
1064 [ + - ][ + - ]: 8640 : << std::setw(9) << grind_perf << " ";
[ + - ]
1065 : :
1066 : : // Augment one-liner status with output indicators
1067 [ + - ][ + + ]: 2880 : if (fielditer() or fieldtime() or fieldrange()) print << 'f';
[ + - ][ + + ]
[ + - ][ + + ]
1068 [ + + ]: 2880 : if (not (m_it % diag)) print << 'd';
1069 [ + - ][ + + ]: 2880 : if (histiter() or histtime() or histrange()) print << 't';
[ + - ][ + + ]
[ + - ][ + + ]
1070 [ + - ][ + + ]: 2880 : if (integiter() or integtime() or integrange()) print << 'i';
[ + - ][ + + ]
[ + - ][ + + ]
1071 [ - + ]: 2880 : if (m_refined) print << 'h';
1072 [ + - ][ + + ]: 2880 : if (lb() and not finished()) print << 'l';
[ + - ][ + + ]
1073 [ + + ][ + - ]: 2880 : if (not benchmark && (not (m_it % rsfreq) || finished())) print << 'c';
[ + - ][ + + ]
1074 [ + + ][ + - ]: 2880 : if (m_deastarted and deactivate()) {
[ + + ]
1075 : : print << "\te:" << m_dea << '/' << m_nchare;
1076 : : }
1077 [ + + ]: 2880 : if (pre) print << "\tp:" << m_pit;
1078 [ + + ]: 2880 : if (mom) print << "\tm:" << m_mit;
1079 : :
1080 : : print << '\n';
1081 : : }
1082 : 18397 : }
1083 : :
1084 : : #include "NoWarning/discretization.def.h"
|
