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-2025 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 "PDFReducer.hpp"
26 : : #include "XystBuildConfig.hpp"
27 : : #include "Box.hpp"
28 : : #include "Transfer.hpp"
29 : : #include "HoleReducer.hpp"
30 : :
31 : : namespace inciter {
32 : :
33 : : static CkReduction::reducerType PDFMerger;
34 : : static CkReduction::reducerType HoleMerger;
35 : : extern ctr::Config g_cfg;
36 : :
37 : : } // inciter::
38 : :
39 : : using inciter::Discretization;
40 : :
41 : 2576 : Discretization::Discretization(
42 : : std::size_t meshid,
43 : : const std::vector< CProxy_Discretization >& disc,
44 : : const CProxy_Transporter& transporter,
45 : : const tk::CProxy_MeshWriter& meshwriter,
46 : : const tk::UnsMesh::CoordMap& coordmap,
47 : : const tk::UnsMesh::Chunk& el,
48 : : const std::map< int, std::unordered_set< std::size_t > >& nodeCommMap,
49 : 2576 : int nc ) :
50 : 2576 : m_meshid( meshid ),
51 : 2576 : m_nchare( nc ),
52 : 2576 : m_it( 0 ),
53 : 2576 : m_itr( 0 ),
54 : 2576 : m_itf( 0 ),
55 : 2576 : m_initial( 1 ),
56 : 2576 : m_t( g_cfg.get< tag::t0 >() ),
57 : 2576 : m_lastDumpTime( -std::numeric_limits< tk::real >::max() ),
58 : 2576 : m_physFieldFloor( 0.0 ),
59 : 2576 : m_physHistFloor( 0.0 ),
60 : 2576 : m_physIntegFloor( 0.0 ),
61 [ + - ]: 2576 : m_rangeFieldFloor( g_cfg.get< tag::fieldout, tag::range >().size(), 0.0 ),
62 [ + - ]: 2576 : m_rangeHistFloor( g_cfg.get< tag::histout, tag::range >().size(), 0.0 ),
63 [ + - ]: 2576 : m_rangeIntegFloor( g_cfg.get< tag::integout, tag::range >().size(), 0.0 ),
64 : 2576 : m_dt( g_cfg.get< tag::dt >() ),
65 : 2576 : m_dtn( m_dt ),
66 : 2576 : m_nvol( 0 ),
67 [ + - ]: 2576 : m_disc( disc ),
68 [ + - ]: 2576 : m_transporter( transporter ),
69 [ + - ]: 2576 : m_meshwriter( meshwriter ),
70 [ + - ]: 2576 : m_el( el ), // fills m_inpoel, m_gid, m_lid
71 : 2576 : m_coord( setCoord( coordmap ) ),
72 : 2576 : m_meshvol( 0.0 ),
73 [ + - ]: 2576 : m_v( m_gid.size(), 0.0 ),
74 [ + - ]: 2576 : m_vol( m_gid.size(), 0.0 ),
75 : 2576 : m_refined( 0 ),
76 : 2576 : m_prevstatus( std::chrono::high_resolution_clock::now() ),
77 : 2576 : m_nrestart( 0 ),
78 : 2576 : m_res( 0.0 ),
79 : 2576 : m_res0( 0.0 ),
80 : 2576 : m_res1( 0.0 ),
81 : 2576 : m_dea( 0 ),
82 : 2576 : m_deastarted( 0 ),
83 : 2576 : m_pit( 0 ),
84 [ + - ]: 15456 : m_mit( 0 )
85 : : // *****************************************************************************
86 : : // Constructor
87 : : //! \param[in] meshid Mesh ID
88 : : //! \param[in] disc Discretization proxy for all meshes
89 : : //! \param[in] transporter Host (Transporter) proxy
90 : : //! \param[in] meshwriter Mesh writer proxy
91 : : //! \param[in] coordmap Coordinates of mesh nodes and their global IDs
92 : : //! \param[in] el Elements of the mesh chunk we operate on
93 : : //! \param[in] nodeCommMap Node lists associated to chare IDs bordering the
94 : : //! mesh chunk we operate on
95 : : //! \param[in] nc Total number of Discretization chares
96 : : // *****************************************************************************
97 : : {
98 [ - + ][ - - ]: 2576 : Assert( !m_inpoel.empty(), "No elements assigned to Discretization chare" );
[ - - ][ - - ]
99 [ + - ][ - + ]: 2576 : Assert( tk::positiveJacobians( m_inpoel, m_coord ),
[ - - ][ - - ]
[ - - ]
100 : : "Jacobian in input mesh to Discretization non-positive" );
101 : : #if not defined(__INTEL_COMPILER) || defined(NDEBUG)
102 : : // The above ifdef skips running the conformity test with the intel compiler
103 : : // in debug mode only. This is necessary because in tk::conforming(), filling
104 : : // up the map can fail with some meshes (only in parallel), e.g., tube.exo,
105 : : // used by some regression tests, due to the intel compiler generating some
106 : : // garbage incorrect code - only in debug, only in parallel, only with that
107 : : // mesh.
108 [ + - ][ - + ]: 2576 : Assert( tk::conforming( m_inpoel, m_coord ),
[ - - ][ - - ]
[ - - ]
109 : : "Input mesh to Discretization not conforming" );
110 : : #endif
111 : :
112 : : // Store node communication map
113 [ + - ][ + - ]: 27504 : for (const auto& [c,map] : nodeCommMap) m_nodeCommMap[c] = map;
[ + + ]
114 : :
115 : : // Get ready for computing/communicating nodal volumes
116 [ + - ]: 2576 : startvol();
117 : :
118 : : // Find host elements of user-specified points where time histories are
119 : : // saved, and save the shape functions evaluated at the point locations
120 : 2576 : const auto& pt = g_cfg.get< tag::histout, tag::points >();
121 [ + + ]: 2816 : for (std::size_t p=0; p<pt.size(); ++p) {
122 : : std::array< tk::real, 4 > N;
123 : 240 : const auto& l = pt[p];
124 [ + + ]: 117607 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
125 [ + - ][ + + ]: 117442 : if (tk::intet( m_coord, m_inpoel, l, e, N )) {
126 [ + - ][ + - ]: 75 : m_histdata.push_back( HistData{{ "p"+std::to_string(p+1), e, N }} );
[ + - ]
127 : 75 : break;
128 : : }
129 : : }
130 : : }
131 : :
132 : : // Register with mesh-transfer (if coupled)
133 [ + - ]: 2576 : if (m_disc.size() == 1) {
134 [ + - ]: 2576 : transfer_initialized();
135 : : } else {
136 [ - - ]: 0 : if (thisIndex == 0) {
137 [ - - ][ - - ]: 0 : transfer::addMesh( thisProxy, m_nchare,
138 [ - - ][ - - ]: 0 : CkCallback(CkIndex_Discretization::transfer_initialized(), thisProxy) );
139 : : }
140 : : }
141 : 2576 : }
142 : :
143 : : void
144 : 2576 : Discretization::transfer_initialized()
145 : : // *****************************************************************************
146 : : // Our mesh has been registered with the mesh-to-mesh transfer (if coupled)
147 : : // *****************************************************************************
148 : : {
149 : : // Compute number of mesh points owned
150 : 2576 : std::size_t npoin = m_gid.size();
151 [ + - ][ + + ]: 318543 : for (auto g : m_gid) if (tk::slave( m_nodeCommMap, g, thisIndex ) ) --npoin;
[ + + ]
152 : :
153 : : // Tell the RTS that the Discretization chares have been created and compute
154 : : // the total number of mesh points across the distributed mesh
155 [ + - ]: 2576 : std::vector< std::size_t > meshdata{ m_meshid, npoin };
156 [ + - ]: 2576 : contribute( meshdata, CkReduction::sum_ulong,
157 [ + - ][ + - ]: 5152 : CkCallback( CkReductionTarget(Transporter,disccreated), m_transporter ) );
158 : 2576 : }
159 : :
160 : : void
161 : 9755 : Discretization::transfer( tk::Fields& u, CkCallback c, bool trflag )
162 : : // *****************************************************************************
163 : : // Initiate solution transfer from background to overset mesh (if coupled)
164 : : // *****************************************************************************
165 : : {
166 [ + - ]: 9755 : if (m_disc.size() == 1) { // not coupled
167 : :
168 : 9755 : c.send();
169 : :
170 : : }
171 : : else {
172 : :
173 : 0 : m_transfer_complete = c;
174 : 0 : m_transfer_sol = static_cast< tk::Fields* >( &u );
175 : 0 : m_trflag = trflag;
176 : :
177 : : // Initiate transfer in 'to' direction
178 [ - - ]: 0 : if (!m_meshid) {
179 [ - - ]: 0 : transfer::setSourceTets( thisProxy, thisIndex, m_inpoel, m_coord, u,
180 [ - - ]: 0 : m_transfer_flag, /*transfer direction=*/ 0,
181 [ - - ][ - - ]: 0 : CkCallback( CkIndex_Discretization::transfer_from(),
182 [ - - ]: 0 : thisProxy[thisIndex] ) );
183 : : }
184 : : else {
185 [ - - ]: 0 : transfer::setDestPoints( thisProxy, thisIndex, m_coord, u,
186 [ - - ]: 0 : m_transfer_flag, m_trflag, /*transfer direction=*/ 0,
187 [ - - ][ - - ]: 0 : CkCallback( CkIndex_Discretization::transfer_from(),
188 [ - - ]: 0 : thisProxy[thisIndex] ) );
189 : : }
190 : :
191 : : }
192 : 9755 : }
193 : :
194 : : void
195 : 0 : Discretization::transfer_from()
196 : : // *****************************************************************************
197 : : // Initiate solution transfer from overset to background mesh
198 : : // *****************************************************************************
199 : : {
200 [ - - ]: 0 : if (!m_meshid) {
201 [ - - ]: 0 : transfer::setDestPoints( thisProxy, thisIndex, m_coord, *m_transfer_sol,
202 : 0 : m_transfer_flag, m_trflag, /*transfer direction=*/ 1,
203 [ - - ]: 0 : m_transfer_complete );
204 : : }
205 : : else {
206 [ - - ]: 0 : transfer::setSourceTets( thisProxy, thisIndex, m_inpoel, m_coord,
207 : 0 : *m_transfer_sol, m_transfer_flag, /*transfer direction=*/ 1,
208 [ - - ]: 0 : m_transfer_complete );
209 : : }
210 : 0 : }
211 : :
212 : : void
213 : 385 : Discretization::intergrid(
214 : : const std::map< int, std::vector< std::size_t > >& bnode )
215 : : // *****************************************************************************
216 : : // Prepare integrid-boundary data structures (if coupled)
217 : : //! \param[in] bnode Boundary-node lists mapped to side sets used in input file
218 : : // *****************************************************************************
219 : : {
220 : 385 : bool multi = g_cfg.get< tag::input >().size() > 1;
221 [ + - ]: 385 : if (!multi) return;
222 [ - - ]: 0 : m_transfer_flag.resize( m_coord[0].size(), -2.0 );
223 [ - - ]: 0 : if (!m_meshid) return;
224 : :
225 : : // Access intergrid-boundary side set ids for this mesh
226 : 0 : const auto& setids = g_cfg.get< tag::overset, tag::intergrid_ >()[ m_meshid ];
227 : :
228 : : // Compile unique set of intergrid-boundary side set ids for this mesh
229 [ - - ]: 0 : std::unordered_set< std::size_t > ibs( begin(setids), end(setids) );
230 [ - - ]: 0 : if (ibs.empty()) return;
231 : :
232 : : // Flag points on intergrid boundary
233 : 0 : std::unordered_set< std::size_t > bp; // points flagged
234 [ - - ]: 0 : for (const auto& [setid,n] : bnode) {
235 [ - - ][ - - ]: 0 : if (ibs.count( static_cast<std::size_t>(setid) )) {
236 [ - - ]: 0 : for (auto g : n) {
237 [ - - ]: 0 : auto i = tk::cref_find(m_lid,g);
238 : 0 : m_transfer_flag[i] = 2.0;
239 [ - - ]: 0 : bp.insert(i);
240 : : }
241 : : }
242 : : }
243 : :
244 : : // Configure number of layers
245 : 0 : auto eps = std::numeric_limits< tk::real >::epsilon();
246 : 0 : const auto& layers = g_cfg.get< tag::overset, tag::layers_ >()[ m_meshid ];
247 : 0 : std::uint64_t ilays = 2; // default number of inner layers
248 : 0 : std::uint64_t mlays = 4; // default number of middle layers
249 : 0 : std::uint64_t olays = 2; // default number of outer layers
250 [ - - ]: 0 : if (layers.size() == 3) {
251 : 0 : ilays = layers[0];
252 : 0 : mlays = layers[1];
253 : 0 : olays = layers[2];
254 : : }
255 : :
256 : : // Add some layers to intergrid boundary
257 [ - - ][ - - ]: 0 : auto psup = tk::genPsup( m_inpoel, 4, tk::genEsup( m_inpoel, 4 ) );
258 [ - - ]: 0 : for (std::uint64_t n=0; n<ilays; ++n) {
259 : 0 : std::unordered_set< std::size_t > add;
260 [ - - ]: 0 : for (auto p : bp) {
261 [ - - ]: 0 : for (auto q : tk::Around(psup,p)) {
262 [ - - ]: 0 : if (std::abs(m_transfer_flag[q]+2.0) < eps) m_transfer_flag[q] = 1.0;
263 [ - - ]: 0 : add.insert(q);
264 : : }
265 : : }
266 [ - - ]: 0 : bp.merge( add );
267 : 0 : add.clear();
268 : 0 : }
269 : :
270 : : // Mark middle layers for buffering / relaxation
271 [ - - ]: 0 : for (std::uint64_t n=0; n<mlays; ++n) {
272 : 0 : std::unordered_set< std::size_t > add;
273 [ - - ]: 0 : for (auto p : bp) {
274 [ - - ]: 0 : for (auto q : tk::Around(psup,p)) {
275 [ - - ]: 0 : if (std::abs(m_transfer_flag[q]+2.0) < eps) m_transfer_flag[q] = -1.0;
276 [ - - ]: 0 : add.insert(q);
277 : : }
278 : : }
279 [ - - ]: 0 : bp.merge( add );
280 : 0 : add.clear();
281 : 0 : }
282 : :
283 : : // Mark outer layers for solution transfer from background to overset
284 [ - - ]: 0 : for (std::uint64_t n=0; n<olays; ++n) {
285 : 0 : std::unordered_set< std::size_t > add;
286 [ - - ]: 0 : for (auto p : bp) {
287 [ - - ]: 0 : for (auto q : tk::Around(psup,p)) {
288 [ - - ]: 0 : if (std::abs(m_transfer_flag[q]+2.0) < eps) m_transfer_flag[q] = 0.0;
289 [ - - ]: 0 : add.insert(q);
290 : : }
291 : : }
292 [ - - ]: 0 : bp.merge( add );
293 : 0 : add.clear();
294 : 0 : }
295 [ - - ]: 0 : }
296 : :
297 : : void
298 : 385 : Discretization::hole(
299 : : const std::map< int, std::vector< std::size_t > >& bface,
300 : : const std::vector< std::size_t >& triinpoel,
301 : : CkCallback c )
302 : : // *****************************************************************************
303 : : // Prepare integrid-boundary data structures (if coupled)
304 : : //! \param[in] bface Boundary-face lists mapped to side sets used in input file
305 : : //! \param[in] triinpoel Boundary-face connectivity (local ids)
306 : : //! \param[in] c Call to continue with when done
307 : : // *****************************************************************************
308 : : {
309 : 385 : bool multi = g_cfg.get< tag::input >().size() > 1;
310 [ + - ][ + - ]: 385 : if (!multi) { c.send(); return; }
311 : :
312 : 0 : m_holcont = c;
313 : :
314 : 0 : std::unordered_map< std::size_t, std::vector< tk::real > > hol;
315 : :
316 [ - - ]: 0 : if (m_meshid) {
317 : : // Access intergrid-boundary side set ids for this mesh
318 : 0 : const auto& setids = g_cfg.get< tag::overset, tag::intergrid_ >()[m_meshid];
319 : : // Compile unique set of intergrid-boundary side set ids for this mesh
320 [ - - ]: 0 : std::unordered_set< std::size_t > ibs( begin(setids), end(setids) );
321 : :
322 : 0 : const auto& x = m_coord[0];
323 : 0 : const auto& y = m_coord[1];
324 : 0 : const auto& z = m_coord[2];
325 : :
326 : : // Collect faces on intergrid boundary
327 [ - - ]: 0 : auto& h = hol[0]; // a single hole for now (given by multiple side sets)
328 [ - - ]: 0 : for (const auto& [setid,faceids] : bface) {
329 [ - - ][ - - ]: 0 : if (ibs.count( static_cast<std::size_t>(setid) )) {
330 [ - - ]: 0 : for (auto f : faceids) {
331 : 0 : const auto t = triinpoel.data() + f*3;
332 [ - - ]: 0 : h.push_back( x[t[0]] );
333 [ - - ]: 0 : h.push_back( y[t[0]] );
334 [ - - ]: 0 : h.push_back( z[t[0]] );
335 [ - - ]: 0 : h.push_back( x[t[1]] );
336 [ - - ]: 0 : h.push_back( y[t[1]] );
337 [ - - ]: 0 : h.push_back( z[t[1]] );
338 [ - - ]: 0 : h.push_back( x[t[2]] );
339 [ - - ]: 0 : h.push_back( y[t[2]] );
340 [ - - ]: 0 : h.push_back( z[t[2]] );
341 : : }
342 : : }
343 : : }
344 : :
345 : : // overset mesh sends hole-parts to background mesh for assembly
346 : : // (allreduce to all partitions of mesh 0)
347 [ - - ]: 0 : auto stream = serialize( hol );
348 [ - - ]: 0 : contribute( stream.first, stream.second.get(), HoleMerger,
349 [ - - ][ - - ]: 0 : CkCallback(CkIndex_Discretization::aggregateHoles(nullptr), m_disc[0]) );
350 : 0 : }
351 : 0 : }
352 : :
353 : : void
354 : 0 : Discretization::aggregateHoles( CkReductionMsg* msg )
355 : : // *****************************************************************************
356 : : // Receive hole data from other meshes
357 : : //! \param[in] msg Aggregated hole data
358 : : //! \note Only background mesh (mesh 0) is supposed to call this.
359 : : // *****************************************************************************
360 : : {
361 : 0 : std::unordered_map< std::size_t, std::vector< tk::real > > inhol;
362 : :
363 : 0 : PUP::fromMem creator( msg->getData() );
364 [ - - ]: 0 : creator | inhol;
365 [ - - ][ - - ]: 0 : delete msg;
366 : :
367 [ - - ]: 0 : for (auto&& [hid,data] : inhol) {
368 [ - - ]: 0 : auto& hole = m_transfer_hole[ hid ];
369 [ - - ][ - - ]: 0 : std::move( begin(data), end(data), std::back_inserter(hole) );
370 : : }
371 : :
372 : : // back to sender overset mesh so it can continue (enough to this partition)
373 [ - - ][ - - ]: 0 : m_disc[ 1 ][ thisIndex ].holeComplete();
374 : :
375 : : // background mesh also complete
376 [ - - ]: 0 : m_holcont.send();
377 : 0 : }
378 : :
379 : : void
380 : 0 : Discretization::holeComplete()
381 : : // *****************************************************************************
382 : : // Hole communication complete
383 : : // *****************************************************************************
384 : : {
385 : 0 : m_holcont.send();
386 : 0 : }
387 : :
388 : : void
389 : 385 : Discretization::holefind()
390 : : // *****************************************************************************
391 : : // Find mesh nodes within holes
392 : : // *****************************************************************************
393 : : {
394 : 385 : bool multi = g_cfg.get< tag::input >().size() > 1;
395 [ + - ]: 385 : if (!multi) return;
396 [ - - ][ - - ]: 0 : if (m_meshid or m_transfer_hole.empty()) return;
[ - - ]
397 : :
398 : : // account for hole symmetry if specified
399 : 0 : const auto& sym = g_cfg.get< tag::overset, tag::sym_ >()[ 1 ];
400 : :
401 : : // collect face centers and normals for hole surface
402 : 0 : std::vector< tk::real > face;
403 [ - - ]: 0 : for (const auto& [hid,tricoord] : m_transfer_hole) { // for each hole
404 [ - - ]: 0 : for (std::size_t t=0; t<tricoord.size()/9; ++t) { // each hole triangle
405 : 0 : const auto f = tricoord.data() + t*9;
406 : 0 : auto x0 = f[0];
407 : 0 : auto y0 = f[1];
408 : 0 : auto z0 = f[2];
409 : 0 : auto x1 = f[3];
410 : 0 : auto y1 = f[4];
411 : 0 : auto z1 = f[5];
412 : 0 : auto x2 = f[6];
413 : 0 : auto y2 = f[7];
414 : 0 : auto z2 = f[8];
415 : 0 : auto cx = (x0 + x1 + x2) / 3.0;
416 : 0 : auto cy = (y0 + y1 + y2) / 3.0;
417 : 0 : auto cz = (z0 + z1 + z2) / 3.0;
418 : : const std::array< tk::real, 3 >
419 : 0 : ba{ x1-x0, y1-y0, z1-z0 }, ca{ x2-x0, y2-y0, z2-z0 };
420 : 0 : auto [nx,ny,nz] = tk::cross( ba, ca );
421 : 0 : nx /= -2.0;
422 : 0 : ny /= -2.0;
423 : 0 : nz /= -2.0;
424 [ - - ]: 0 : face.push_back( cx );
425 [ - - ]: 0 : face.push_back( cy );
426 [ - - ]: 0 : face.push_back( cz );
427 [ - - ]: 0 : face.push_back( nx );
428 [ - - ]: 0 : face.push_back( ny );
429 [ - - ]: 0 : face.push_back( nz );
430 [ - - ]: 0 : if (sym == "z") { // augment hole with symmetric part
431 [ - - ]: 0 : face.push_back( cx );
432 [ - - ]: 0 : face.push_back( cy );
433 [ - - ]: 0 : face.push_back( -cz );
434 [ - - ]: 0 : face.push_back( nx );
435 [ - - ]: 0 : face.push_back( ny );
436 [ - - ]: 0 : face.push_back( -nz );
437 : : }
438 : : }
439 : : }
440 : :
441 : 0 : const auto& x = m_coord[0];
442 : 0 : const auto& y = m_coord[1];
443 : 0 : const auto& z = m_coord[2];
444 : 0 : auto npoin = x.size();
445 : :
446 : : // compute integral for finding hole nodes on background mesh
447 : 0 : auto inteps = 1.0;
448 : 0 : auto intval = 4.0 * M_PI;
449 : 0 : std::unordered_set< std::size_t > hp; // hole points
450 [ - - ]: 0 : for (std::size_t i=0; i<npoin; ++i) {
451 : 0 : tk::real holeint = 0.0;
452 [ - - ]: 0 : for (std::size_t t=0; t<face.size()/6; ++t) {
453 : 0 : const auto f = face.data() + t*6;
454 : 0 : auto dx = f[0] - x[i];
455 : 0 : auto dy = f[1] - y[i];
456 : 0 : auto dz = f[2] - z[i];
457 : 0 : auto r = std::pow( dx*dx + dy*dy + dz*dz, 1.5 );
458 : 0 : auto vx = dx / r;
459 : 0 : auto vy = dy / r;
460 : 0 : auto vz = dz / r;
461 : 0 : holeint += vx*f[3] + vy*f[4] + vz*f[5];
462 : : }
463 [ - - ]: 0 : if (std::abs(holeint - intval) < inteps) {
464 : 0 : m_transfer_flag[i] = 2.0;
465 [ - - ]: 0 : hp.insert(i);
466 : : }
467 : : }
468 : :
469 : : // grow hole by one extra layer
470 [ - - ][ - - ]: 0 : auto psup = tk::genPsup( m_inpoel, 4, tk::genEsup( m_inpoel, 4 ) );
471 : 0 : auto eps = std::numeric_limits< tk::real >::epsilon();
472 [ - - ]: 0 : for (auto p : hp) {
473 [ - - ]: 0 : for (auto q : tk::Around(psup,p)) {
474 [ - - ]: 0 : if (std::abs(m_transfer_flag[q]+2.0) < eps) m_transfer_flag[q] = 2.0;
475 : : }
476 : : }
477 : 0 : }
478 : :
479 : : void
480 : 0 : Discretization::resizePostAMR(
481 : : const tk::UnsMesh::Chunk& chunk,
482 : : const tk::UnsMesh::Coords& coord,
483 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
484 : : nodeCommMap,
485 : : const std::set< std::size_t >& /*removedNodes*/ )
486 : : // *****************************************************************************
487 : : // Resize mesh data structures after mesh refinement
488 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
489 : : //! \param[in] coord New mesh node coordinates
490 : : //! \param[in] nodeCommMap New node communication map
491 : : //! \param[in] removedNodes Newly removed mesh node local ids
492 : : // *****************************************************************************
493 : : {
494 : 0 : m_el = chunk; // updates m_inpoel, m_gid, m_lid
495 : 0 : m_nodeCommMap = nodeCommMap;
496 : :
497 : : // Update mesh volume container size
498 [ - - ]: 0 : m_vol.resize( m_gid.size(), 0.0 );
499 : :
500 : : // update mesh node coordinates
501 : 0 : m_coord = coord;
502 : 0 : }
503 : :
504 : : void
505 : 2576 : Discretization::startvol()
506 : : // *****************************************************************************
507 : : // Get ready for (re-)computing/communicating nodal volumes
508 : : // *****************************************************************************
509 : : {
510 : 2576 : m_nvol = 0;
511 [ + - ][ + - ]: 2576 : thisProxy[ thisIndex ].wait4vol();
512 : :
513 : : // Zero out mesh volume container
514 [ + - ]: 2576 : std::fill( begin(m_vol), end(m_vol), 0.0 );
515 : :
516 : : // Clear receive buffer that will be used for collecting nodal volumes
517 : 2576 : m_volc.clear();
518 : 2576 : }
519 : :
520 : : void
521 : 804 : Discretization::registerReducers()
522 : : // *****************************************************************************
523 : : // Configure Charm++ reduction types
524 : : //! \details Since this is a [initnode] routine, see the .ci file, the
525 : : //! Charm++ runtime system executes the routine exactly once on every
526 : : //! logical node early on in the Charm++ init sequence. Must be static as
527 : : //! it is called without an object. See also: Section "Initializations at
528 : : //! Program Startup" at in the Charm++ manual
529 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
530 : : // *****************************************************************************
531 : : {
532 : 804 : PDFMerger = CkReduction::addReducer( tk::mergeUniPDFs );
533 : 804 : HoleMerger = CkReduction::addReducer( mergeHole );
534 : 804 : }
535 : :
536 : : tk::UnsMesh::Coords
537 : 2576 : Discretization::setCoord( const tk::UnsMesh::CoordMap& coordmap )
538 : : // *****************************************************************************
539 : : // Set mesh coordinates based on coordinates map
540 : : // *****************************************************************************
541 : : {
542 [ - + ][ - - ]: 2576 : Assert( coordmap.size() == m_gid.size(), "Size mismatch" );
[ - - ][ - - ]
543 [ - + ][ - - ]: 2576 : Assert( coordmap.size() == m_lid.size(), "Size mismatch" );
[ - - ][ - - ]
544 : :
545 : 2576 : tk::UnsMesh::Coords coord;
546 [ + - ]: 2576 : coord[0].resize( coordmap.size() );
547 [ + - ]: 2576 : coord[1].resize( coordmap.size() );
548 [ + - ]: 2576 : coord[2].resize( coordmap.size() );
549 : :
550 [ + + ]: 318543 : for (const auto& [ gid, coords ] : coordmap) {
551 [ + - ]: 315967 : auto i = tk::cref_find( m_lid, gid );
552 : 315967 : coord[0][i] = coords[0];
553 : 315967 : coord[1][i] = coords[1];
554 : 315967 : coord[2][i] = coords[2];
555 : : }
556 : :
557 : 2576 : return coord;
558 : 0 : }
559 : :
560 : : void
561 : 1682 : Discretization::remap(
562 : : const std::unordered_map< std::size_t, std::size_t >& map )
563 : : // *****************************************************************************
564 : : // Remap mesh data based on new local ids
565 : : //! \param[in] map Mapping of old->new local ids
566 : : // *****************************************************************************
567 : : {
568 : : // Remap connectivity containing local IDs
569 [ + - ][ + + ]: 2787754 : for (auto& l : m_inpoel) l = tk::cref_find(map,l);
570 : :
571 : : // Remap global->local id map
572 [ + - ][ + + ]: 200013 : for (auto& [g,l] : m_lid) l = tk::cref_find(map,l);
573 : :
574 : : // Remap global->local id map
575 : 1682 : auto maxid = std::numeric_limits< std::size_t >::max();
576 [ + - ]: 1682 : std::vector< std::size_t > newgid( m_gid.size(), maxid );
577 [ + + ]: 200013 : for (const auto& [o,n] : map) newgid[n] = m_gid[o];
578 : 1682 : m_gid = std::move( newgid );
579 : :
580 [ + - ][ - + ]: 200013 : Assert( std::all_of( m_gid.cbegin(), m_gid.cend(),
[ - - ][ - - ]
[ - - ]
581 : : [=](std::size_t i){ return i < maxid; } ),
582 : : "Not all gid have been remapped" );
583 : :
584 : : // Remap nodal volumes (with contributions along chare-boundaries)
585 [ + - ]: 1682 : std::vector< tk::real > newvol( m_vol.size(), 0.0 );
586 [ + + ]: 200013 : for (const auto& [o,n] : map) newvol[n] = m_vol[o];
587 : 1682 : m_vol = std::move( newvol );
588 : :
589 : : // Remap nodal volumes (without contributions along chare-boundaries)
590 [ + - ]: 1682 : std::vector< tk::real > newv( m_v.size(), 0.0 );
591 [ + + ]: 200013 : for (const auto& [o,n] : map) newv[n] = m_v[o];
592 : 1682 : m_v = std::move( newv );
593 : :
594 : : // Remap locations of node coordinates
595 : 1682 : tk::UnsMesh::Coords newcoord;
596 : 1682 : auto npoin = m_coord[0].size();
597 [ + - ]: 1682 : newcoord[0].resize( npoin );
598 [ + - ]: 1682 : newcoord[1].resize( npoin );
599 [ + - ]: 1682 : newcoord[2].resize( npoin );
600 [ + + ]: 200013 : for (const auto& [o,n] : map) {
601 : 198331 : newcoord[0][n] = m_coord[0][o];
602 : 198331 : newcoord[1][n] = m_coord[1][o];
603 : 198331 : newcoord[2][n] = m_coord[2][o];
604 : : }
605 : 1682 : m_coord = std::move( newcoord );
606 : 1682 : }
607 : :
608 : : void
609 : 2576 : Discretization::setRefiner( const CProxy_Refiner& ref )
610 : : // *****************************************************************************
611 : : // Set Refiner Charm++ proxy
612 : : //! \param[in] ref Incoming refiner proxy to store
613 : : // *****************************************************************************
614 : : {
615 : 2576 : m_refiner = ref;
616 : 2576 : }
617 : :
618 : : void
619 : 2576 : Discretization::vol()
620 : : // *****************************************************************************
621 : : // Sum mesh volumes to nodes, start communicating them on chare-boundaries
622 : : // *****************************************************************************
623 : : {
624 : 2576 : const auto& x = m_coord[0];
625 : 2576 : const auto& y = m_coord[1];
626 : 2576 : const auto& z = m_coord[2];
627 : :
628 : : // Compute nodal volumes on our chunk of the mesh
629 [ + + ]: 1138146 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
630 : 1135570 : const auto N = m_inpoel.data() + e*4;
631 : : const std::array< tk::real, 3 >
632 : 1135570 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
633 : 1135570 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
634 : 1135570 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
635 : 1135570 : const auto J = tk::triple( ba, ca, da ) / 24.0;
636 [ - + ][ - - ]: 1135570 : ErrChk( J > 0, "Element Jacobian non-positive: PE:" +
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
637 : : std::to_string(CkMyPe()) + ", node IDs: " +
638 : : std::to_string(m_gid[N[0]]) + ',' +
639 : : std::to_string(m_gid[N[1]]) + ',' +
640 : : std::to_string(m_gid[N[2]]) + ',' +
641 : : std::to_string(m_gid[N[3]]) + ", coords: (" +
642 : : std::to_string(x[N[0]]) + ", " +
643 : : std::to_string(y[N[0]]) + ", " +
644 : : std::to_string(z[N[0]]) + "), (" +
645 : : std::to_string(x[N[1]]) + ", " +
646 : : std::to_string(y[N[1]]) + ", " +
647 : : std::to_string(z[N[1]]) + "), (" +
648 : : std::to_string(x[N[2]]) + ", " +
649 : : std::to_string(y[N[2]]) + ", " +
650 : : std::to_string(z[N[2]]) + "), (" +
651 : : std::to_string(x[N[3]]) + ", " +
652 : : std::to_string(y[N[3]]) + ", " +
653 : : std::to_string(z[N[3]]) + ')' );
654 : : // scatter add V/4 to nodes
655 [ + + ]: 5677850 : for (std::size_t j=0; j<4; ++j) m_vol[N[j]] += J;
656 : : }
657 : :
658 : : // Store nodal volumes without contributions from other chares on
659 : : // chare-boundaries
660 : 2576 : m_v = m_vol;
661 : :
662 : : // Send our nodal volume contributions to neighbor chares
663 [ + + ]: 2576 : if (m_nodeCommMap.empty()) {
664 : 78 : comvol_complete();
665 : : } else {
666 [ + + ]: 27426 : for (const auto& [c,n] : m_nodeCommMap) {
667 [ + - ]: 24928 : std::vector< tk::real > v( n.size() );
668 : 24928 : std::size_t j = 0;
669 [ + - ][ + + ]: 176660 : for (auto i : n) v[ j++ ] = m_vol[ tk::cref_find(m_lid,i) ];
670 [ + - ][ + - ]: 49856 : thisProxy[c].comvol( thisIndex,
671 [ + - ]: 49856 : std::vector<std::size_t>(begin(n), end(n)), v );
672 : 24928 : }
673 : : }
674 : :
675 : 2576 : ownvol_complete();
676 : 2576 : }
677 : :
678 : : void
679 : 24928 : Discretization::comvol( int c,
680 : : const std::vector< std::size_t >& gid,
681 : : const std::vector< tk::real >& nodevol )
682 : : // *****************************************************************************
683 : : // Receive nodal volumes on chare-boundaries
684 : : //! \param[in] c Sender chare id
685 : : //! \param[in] gid Global mesh node IDs at which we receive volume contributions
686 : : //! \param[in] nodevol Partial sums of nodal volume contributions to
687 : : //! chare-boundary nodes
688 : : // *****************************************************************************
689 : : {
690 [ - + ][ - - ]: 24928 : Assert( nodevol.size() == gid.size(), "Size mismatch" );
[ - - ][ - - ]
691 : :
692 : 24928 : auto& cvolc = m_cvolc[c];
693 [ + + ]: 176660 : for (std::size_t i=0; i<gid.size(); ++i) {
694 : 151732 : m_volc[ gid[i] ] += nodevol[i];
695 : 151732 : cvolc[ tk::cref_find(m_lid,gid[i]) ] = nodevol[i];
696 : : }
697 : :
698 [ + + ]: 24928 : if (++m_nvol == m_nodeCommMap.size()) {
699 : 2498 : m_nvol = 0;
700 : 2498 : comvol_complete();
701 : : }
702 : 24928 : }
703 : :
704 : : void
705 : 2576 : Discretization::totalvol()
706 : : // *****************************************************************************
707 : : // Sum mesh volumes and contribute own mesh volume to total volume
708 : : // *****************************************************************************
709 : : {
710 : : // Add received contributions to nodal volumes
711 [ + + ]: 93337 : for (const auto& [gid, vol] : m_volc)
712 [ + - ]: 90761 : m_vol[ tk::cref_find(m_lid,gid) ] += vol;
713 : :
714 : : // Clear receive buffer
715 : 2576 : tk::destroy(m_volc);
716 : :
717 : : // Sum mesh volume to host
718 : : std::vector< tk::real > tvol{ 0.0,
719 : 2576 : static_cast<tk::real>(m_initial),
720 [ + - ]: 2576 : static_cast<tk::real>(m_meshid) };
721 [ + + ]: 318543 : for (auto v : m_v) tvol[0] += v;
722 [ + - ]: 2576 : contribute( tvol, CkReduction::sum_double,
723 [ + - ][ + - ]: 5152 : CkCallback(CkReductionTarget(Transporter,totalvol), m_transporter) );
724 : 2576 : }
725 : :
726 : : void
727 : 2576 : Discretization::stat( tk::real mesh_volume )
728 : : // *****************************************************************************
729 : : // Compute mesh cell statistics
730 : : //! \param[in] mesh_volume Total mesh volume
731 : : // *****************************************************************************
732 : : {
733 : : // Store total mesh volume
734 : 2576 : m_meshvol = mesh_volume;
735 : :
736 : 2576 : const auto& x = m_coord[0];
737 : 2576 : const auto& y = m_coord[1];
738 : 2576 : const auto& z = m_coord[2];
739 : :
740 : 2576 : auto MIN = -std::numeric_limits< tk::real >::max();
741 : 2576 : auto MAX = std::numeric_limits< tk::real >::max();
742 [ + - ]: 2576 : std::vector< tk::real > min( 6, MAX );
743 [ + - ]: 2576 : std::vector< tk::real > max( 6, MIN );
744 [ + - ]: 2576 : std::vector< tk::real > sum( 9, 0.0 );
745 : 2576 : tk::UniPDF edgePDF( 1e-4 );
746 : 2576 : tk::UniPDF volPDF( 1e-4 );
747 : 2576 : tk::UniPDF ntetPDF( 1e-4 );
748 : :
749 : : // Compute points surrounding points
750 [ + - ][ + - ]: 2576 : auto psup = tk::genPsup( m_inpoel, 4, tk::genEsup(m_inpoel,4) );
751 [ - + ][ - - ]: 2576 : Assert( psup.second.size()-1 == m_gid.size(),
[ - - ][ - - ]
752 : : "Number of mesh points and number of global IDs unequal" );
753 : :
754 : : // Compute edge length statistics
755 : : // Note that while the min and max edge lengths are independent of the number
756 : : // of CPUs (by the time they are aggregated across all chares), the sum of
757 : : // the edge lengths and the edge length PDF are not. This is because the
758 : : // edges on the chare-boundary are counted multiple times and we
759 : : // conscientiously do not make an effort to precisely compute this, because
760 : : // that would require communication and more complex logic. Since these
761 : : // statistics are intended as simple average diagnostics, we ignore these
762 : : // small differences. For reproducible average edge lengths and edge length
763 : : // PDFs, run the mesh in serial.
764 : 2576 : tk::UnsMesh::EdgeSet edges;
765 [ + + ]: 318543 : for (std::size_t p=0; p<m_gid.size(); ++p) {
766 [ + + ]: 3620953 : for (auto i : tk::Around(psup,p)) {
767 : 3304986 : const auto dx = x[ i ] - x[ p ];
768 : 3304986 : const auto dy = y[ i ] - y[ p ];
769 : 3304986 : const auto dz = z[ i ] - z[ p ];
770 : 3304986 : const auto length = std::sqrt( dx*dx + dy*dy + dz*dz );
771 [ + + ]: 3304986 : if (length < min[0]) min[0] = length;
772 [ + + ]: 3304986 : if (length > max[0]) max[0] = length;
773 : 3304986 : sum[0] += 1.0;
774 : 3304986 : sum[1] += length;
775 [ + - ]: 3304986 : edgePDF.add( length );
776 [ + - ]: 3304986 : edges.insert( { m_gid[i], m_gid[p] } );
777 : : }
778 : : }
779 : :
780 : : // Compute mesh cell volume statistics
781 [ + + ]: 1138146 : for (std::size_t e=0; e<m_inpoel.size()/4; ++e) {
782 : 1135570 : const std::array< std::size_t, 4 > N{{ m_inpoel[e*4+0], m_inpoel[e*4+1],
783 : 1135570 : m_inpoel[e*4+2], m_inpoel[e*4+3] }};
784 : : const std::array< tk::real, 3 >
785 : 1135570 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
786 : 1135570 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
787 : 1135570 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
788 : 1135570 : const auto L = std::cbrt( tk::triple( ba, ca, da ) / 6.0 );
789 [ + + ]: 1135570 : if (L < min[1]) min[1] = L;
790 [ + + ]: 1135570 : if (L > max[1]) max[1] = L;
791 : 1135570 : sum[2] += 1.0;
792 : 1135570 : sum[3] += L;
793 [ + - ]: 1135570 : volPDF.add( L );
794 : : }
795 : :
796 : : // Contribute statistics
797 : 2576 : sum[4] = 1.0;
798 : 2576 : min[2] = max[2] = sum[5] = static_cast< tk::real >( m_inpoel.size() / 4 );
799 : 2576 : min[3] = max[3] = sum[6] = static_cast< tk::real >( m_gid.size() );
800 : 2576 : min[4] = max[4] = sum[7] = static_cast< tk::real >( edges.size() );
801 : 2576 : min[5] = max[5] = sum[8] =
802 : 2576 : static_cast< tk::real >( tk::sumvalsize(m_nodeCommMap) ) /
803 : 2576 : static_cast< tk::real >( m_gid.size() );
804 [ + - ]: 2576 : ntetPDF.add( min[2] );
805 : :
806 [ + - ]: 2576 : min.push_back( static_cast<tk::real>(m_meshid) );
807 [ + - ]: 2576 : max.push_back( static_cast<tk::real>(m_meshid) );
808 [ + - ]: 2576 : sum.push_back( static_cast<tk::real>(m_meshid) );
809 : :
810 : : // Contribute to mesh statistics across all Discretization chares
811 [ + - ]: 2576 : contribute( min, CkReduction::min_double,
812 [ + - ][ + - ]: 5152 : CkCallback(CkReductionTarget(Transporter,minstat), m_transporter) );
813 [ + - ]: 2576 : contribute( max, CkReduction::max_double,
814 [ + - ][ + - ]: 5152 : CkCallback(CkReductionTarget(Transporter,maxstat), m_transporter) );
815 [ + - ]: 2576 : contribute( sum, CkReduction::sum_double,
816 [ + - ][ + - ]: 5152 : CkCallback(CkReductionTarget(Transporter,sumstat), m_transporter) );
817 : :
818 : : // Serialize PDFs to raw stream
819 [ + - ][ + - ]: 10304 : auto stream = tk::serialize( m_meshid, { edgePDF, volPDF, ntetPDF } );
[ + + ][ - - ]
820 : : // Create Charm++ callback function for reduction of PDFs with
821 : : // Transporter::pdfstat() as the final target where the results will appear.
822 [ + - ][ + - ]: 2576 : CkCallback cb( CkIndex_Transporter::pdfstat(nullptr), m_transporter );
823 : : // Contribute serialized PDF of partial sums to host via Charm++ reduction
824 [ + - ]: 2576 : contribute( stream.first, stream.second.get(), PDFMerger, cb );
825 [ + - ][ + - ]: 5152 : }
[ + - ][ - - ]
[ - - ]
826 : :
827 : : void
828 : 2576 : Discretization::boxvol()
829 : : // *****************************************************************************
830 : : // Compute total box IC volume
831 : : // *****************************************************************************
832 : : {
833 : : // Determine which nodes reside in user-defined IC box(es) if any
834 [ + - ]: 2576 : m_boxnodes = problems::boxnodes( m_coord );
835 : :
836 : : // Compute partial box IC volume (just add up all boxes)
837 : 2576 : tk::real boxvol = 0.0;
838 : : // cppcheck-suppress useStlAlgorithm
839 [ + + ][ + + ]: 4055 : for (const auto& b : m_boxnodes) for (auto i : b) boxvol += m_v[i];
840 : :
841 : : // Sum up box IC volume across all chares
842 [ + - ]: 2576 : std::vector< tk::real > meshdata{ boxvol, static_cast<tk::real>(m_meshid) };
843 [ + - ]: 2576 : contribute( meshdata, CkReduction::sum_double,
844 [ + - ][ + - ]: 5152 : CkCallback(CkReductionTarget(Transporter,boxvol), m_transporter) );
845 : 2576 : }
846 : :
847 : : void
848 : 3547 : Discretization::write(
849 : : const std::vector< std::size_t >& inpoel,
850 : : const tk::UnsMesh::Coords& coord,
851 : : const std::map< int, std::vector< std::size_t > >& bface,
852 : : const std::map< int, std::vector< std::size_t > >& bnode,
853 : : const std::vector< std::size_t >& triinpoel,
854 : : const std::vector< std::string>& elemfieldnames,
855 : : const std::vector< std::string>& nodefieldnames,
856 : : const std::vector< std::string>& elemsurfnames,
857 : : const std::vector< std::string>& nodesurfnames,
858 : : const std::vector< std::vector< tk::real > >& elemfields,
859 : : const std::vector< std::vector< tk::real > >& nodefields,
860 : : const std::vector< std::vector< tk::real > >& elemsurfs,
861 : : const std::vector< std::vector< tk::real > >& nodesurfs,
862 : : CkCallback c )
863 : : // *****************************************************************************
864 : : // Output mesh and fields data (solution dump) to file(s)
865 : : //! \param[in] inpoel Mesh connectivity for the mesh chunk to be written
866 : : //! \param[in] coord Node coordinates of the mesh chunk to be written
867 : : //! \param[in] bface Map of boundary-face lists mapped to corresponding side set
868 : : //! ids for this mesh chunk
869 : : //! \param[in] bnode Map of boundary-node lists mapped to corresponding side set
870 : : //! ids for this mesh chunk
871 : : //! \param[in] triinpoel Interconnectivity of points and boundary-face in this
872 : : //! mesh chunk
873 : : //! \param[in] elemfieldnames Names of element fields to be output to file
874 : : //! \param[in] nodefieldnames Names of node fields to be output to file
875 : : //! \param[in] elemsurfnames Names of elem surface fields to be output to file
876 : : //! \param[in] nodesurfnames Names of node surface fields to be output to file
877 : : //! \param[in] elemfields Field data in mesh elements to output to file
878 : : //! \param[in] nodefields Field data in mesh nodes to output to file
879 : : //! \param[in] elemsurfs Surface field data in mesh elements to output to file
880 : : //! \param[in] nodesurfs Surface field data in mesh nodes to output to file
881 : : //! \param[in] c Function to continue with after the write
882 : : //! \details Since m_meshwriter is a Charm++ chare group, it never migrates and
883 : : //! an instance is guaranteed on every PE. We index the first PE on every
884 : : //! logical compute node. In Charm++'s non-SMP mode, a node is the same as a
885 : : //! PE, so the index is the same as CkMyPe(). In SMP mode the index is the
886 : : //! first PE on every logical node. In non-SMP mode this yields one or more
887 : : //! output files per PE with zero or non-zero virtualization, respectively. If
888 : : //! there are multiple chares on a PE, the writes are serialized per PE, since
889 : : //! only a single entry method call can be executed at any given time. In SMP
890 : : //! mode, still the same number of files are output (one per chare), but the
891 : : //! output is serialized through the first PE of each compute node. In SMP
892 : : //! mode, channeling multiple files via a single PE on each node is required
893 : : //! by NetCDF and HDF5, as well as ExodusII, since none of these libraries are
894 : : //! thread-safe.
895 : : // *****************************************************************************
896 : : {
897 : : // If the previous iteration refined (or moved) the mesh or this is called
898 : : // before the first time step, we also output the mesh.
899 : 3547 : bool meshoutput = m_itf == 0 ? true : false;
900 : :
901 : 3547 : auto eps = std::numeric_limits< tk::real >::epsilon();
902 : 3547 : bool fieldoutput = false;
903 : :
904 : : // Output field data only if there is no dump at this physical time yet
905 [ + - ]: 3547 : if (std::abs(m_lastDumpTime - m_t) > eps ) {
906 : 3547 : m_lastDumpTime = m_t;
907 : 3547 : ++m_itf;
908 : 3547 : fieldoutput = true;
909 : : }
910 : :
911 : 3547 : bool multi = g_cfg.get< tag::input >().size() > 1;
912 [ - + ]: 3547 : const auto& ft = multi ? g_cfg.get< tag::fieldout_ >()[ m_meshid ]
913 : 3547 : : g_cfg.get< tag::fieldout >();
914 : :
915 : 3547 : const auto& f = ft.get< tag::sidesets >();
916 [ + - ]: 3547 : std::set< int > outsets( begin(f), end(f) );
917 : :
918 : : m_meshwriter[ CkNodeFirst( CkMyNode() ) ].
919 [ + - ][ + - ]: 7094 : write( m_meshid, meshoutput, fieldoutput, m_itr, m_itf, m_t, thisIndex,
920 : 3547 : g_cfg.get< tag::output >(),
921 : : inpoel, coord, bface, bnode, triinpoel, elemfieldnames,
922 : : nodefieldnames, elemsurfnames, nodesurfnames, elemfields, nodefields,
923 : : elemsurfs, nodesurfs, outsets, c );
924 : 3547 : }
925 : :
926 : : void
927 : 40576 : Discretization::setdt( tk::real newdt )
928 : : // *****************************************************************************
929 : : // Set time step size
930 : : //! \param[in] newdt Size of the new time step
931 : : // *****************************************************************************
932 : : {
933 : 40576 : m_dtn = m_dt;
934 : 40576 : m_dt = newdt;
935 : :
936 : : // Truncate the size of last time step
937 : 40576 : const auto term = g_cfg.get< tag::term >();
938 [ + + ]: 40576 : if (m_t+m_dt > term) m_dt = term - m_t;
939 : 40576 : }
940 : :
941 : : void
942 : 40608 : Discretization::next()
943 : : // *****************************************************************************
944 : : // Prepare for next step
945 : : // *****************************************************************************
946 : : {
947 : : // Update floor of physics time divided by output interval times
948 : 40608 : const auto eps = std::numeric_limits< tk::real >::epsilon();
949 : 40608 : bool multi = g_cfg.get< tag::input >().size() > 1;
950 [ - + ]: 40608 : const auto& ftab = multi ? g_cfg.get< tag::fieldout_ >()[ m_meshid ]
951 : 40608 : : g_cfg.get< tag::fieldout >();
952 : 40608 : const auto ft = ftab.get< tag::time >();
953 [ + + ]: 40608 : if (ft > eps) m_physFieldFloor = std::floor( m_t / ft );
954 : 40608 : const auto ht = g_cfg.get< tag::histout, tag::time >();
955 [ + - ]: 40608 : if (ht > eps) m_physHistFloor = std::floor( m_t / ht );
956 : 40608 : const auto it = g_cfg.get< tag::integout, tag::time >();
957 [ + - ]: 40608 : if (it > eps) m_physIntegFloor = std::floor( m_t / it );
958 : :
959 : : // Update floors of physics time divided by output interval times for ranges
960 : 40608 : const auto& rf = ftab.get< tag::range >();
961 [ + + ]: 41008 : for (std::size_t i=0; i<rf.size(); ++i) {
962 [ + + ][ + + ]: 400 : if (m_t > rf[i][0] and m_t < rf[i][1]) {
[ + + ]
963 : 32 : m_rangeFieldFloor[i] = std::floor( m_t / rf[i][2] );
964 : : }
965 : : }
966 : :
967 : 40608 : const auto& rh = g_cfg.get< tag::histout, tag::range >();
968 [ + + ]: 42848 : for (std::size_t i=0; i<rh.size(); ++i) {
969 [ + + ][ + + ]: 2240 : if (m_t > rh[i][0] and m_t < rh[i][1]) {
[ + + ]
970 : 620 : m_rangeHistFloor[i] = std::floor( m_t / rh[i][2] );
971 : : }
972 : : }
973 : :
974 : 40608 : const auto& ri = g_cfg.get< tag::integout, tag::range >();
975 [ + + ]: 40648 : for (std::size_t i=0; i<ri.size(); ++i) {
976 [ + + ][ + + ]: 40 : if (m_t > ri[i][0] and m_t < ri[i][1]) {
[ + + ]
977 : 10 : m_rangeIntegFloor[i] = std::floor( m_t / ri[i][2] );
978 : : }
979 : : }
980 : :
981 : 40608 : ++m_it;
982 : 40608 : m_t += m_dt;
983 : 40608 : }
984 : :
985 : : void
986 : 2574 : Discretization::grindZero()
987 : : // *****************************************************************************
988 : : // Zero grind-time
989 : : // *****************************************************************************
990 : : {
991 : 2574 : m_prevstatus = std::chrono::high_resolution_clock::now();
992 : :
993 [ + + ][ + - ]: 2574 : if (thisIndex == 0 && m_meshid == 0) {
994 [ + - ]: 260 : tk::Print() << "Starting time stepping ...\n";
995 : : }
996 : 2574 : }
997 : :
998 : : bool
999 : 38032 : Discretization::restarted( int nrestart )
1000 : : // *****************************************************************************
1001 : : // Detect if just returned from a checkpoint and if so, zero timers
1002 : : //! \param[in] nrestart Number of times restarted
1003 : : //! \return True if restart detected
1004 : : // *****************************************************************************
1005 : : {
1006 : : // Detect if just restarted from checkpoint:
1007 : : // nrestart == -1 if there was no checkpoint this step
1008 : : // m_nrestart == nrestart if there was a checkpoint this step
1009 : : // if both false, just restarted from a checkpoint
1010 [ + + ][ + - ]: 38032 : bool restarted = nrestart != -1 and m_nrestart != nrestart;
1011 : :
1012 : : // If just restarted from checkpoint
1013 [ + + ]: 38032 : if (restarted) {
1014 : : // Update number of restarts
1015 : 30 : m_nrestart = nrestart;
1016 : : // Start timer measuring time stepping wall clock time
1017 : 30 : m_timer.zero();
1018 : : // Zero grind-timer
1019 : 30 : grindZero();
1020 : : }
1021 : :
1022 : 38032 : return restarted;
1023 : : }
1024 : :
1025 : : std::string
1026 : 726 : Discretization::histfilename( const std::string& id, std::streamsize precision )
1027 : : // *****************************************************************************
1028 : : // Construct history output filename
1029 : : //! \param[in] id History point id
1030 : : //! \param[in] precision Floating point precision to use for output
1031 : : //! \return History file name
1032 : : // *****************************************************************************
1033 : : {
1034 [ + - ]: 726 : auto of = g_cfg.get< tag::output >();
1035 [ + - ]: 726 : std::stringstream ss;
1036 : :
1037 [ + - ][ + - ]: 726 : ss << std::setprecision(static_cast<int>(precision)) << of << ".hist." << id;
[ + - ]
1038 : :
1039 [ + - ]: 1452 : return ss.str();
1040 : 726 : }
1041 : :
1042 : : void
1043 : 74 : Discretization::histheader( std::vector< std::string >&& names )
1044 : : // *****************************************************************************
1045 : : // Output headers for time history files (one for each point)
1046 : : //! \param[in] names History output variable names
1047 : : // *****************************************************************************
1048 : : {
1049 [ + + ]: 149 : for (const auto& h : m_histdata) {
1050 : 75 : auto prec = g_cfg.get< tag::histout, tag:: precision >();
1051 [ + - ]: 150 : tk::DiagWriter hw( histfilename( h.get< tag::id >(), prec ),
1052 : 75 : g_cfg.get< tag::histout, tag::format >(),
1053 [ + - ]: 75 : prec );
1054 [ + - ]: 75 : hw.header( names );
1055 : 75 : }
1056 : 74 : }
1057 : :
1058 : : void
1059 : 1238 : Discretization::history( std::vector< std::vector< tk::real > >&& data )
1060 : : // *****************************************************************************
1061 : : // Output time history for a time step
1062 : : //! \param[in] data Time history data for all variables and equations integrated
1063 : : // *****************************************************************************
1064 : : {
1065 [ - + ][ - - ]: 1238 : Assert( data.size() == m_histdata.size(), "Size mismatch" );
[ - - ][ - - ]
1066 : :
1067 : 1238 : std::size_t i = 0;
1068 [ + + ]: 1889 : for (const auto& h : m_histdata) {
1069 : 651 : auto prec = g_cfg.get< tag::histout, tag::precision >();
1070 [ + - ]: 1302 : tk::DiagWriter hw( histfilename( h.get< tag::id >(), prec ),
1071 : 651 : g_cfg.get< tag::histout, tag::format >(),
1072 : : prec,
1073 [ + - ]: 651 : std::ios_base::app );
1074 [ + - ]: 651 : hw.write( m_it, m_t, m_dt, data[i] );
1075 : 651 : ++i;
1076 : 651 : }
1077 : 1238 : }
1078 : :
1079 : : bool
1080 : 43668 : Discretization::fielditer() const
1081 : : // *****************************************************************************
1082 : : // Decide if field output iteration count interval is hit
1083 : : //! \return True if field output iteration count interval is hit
1084 : : // *****************************************************************************
1085 : : {
1086 [ + + ]: 43668 : if (g_cfg.get< tag::benchmark >()) return false;
1087 : :
1088 : 25068 : bool multi = g_cfg.get< tag::input >().size() > 1;
1089 [ - + ]: 25068 : const auto& ft = multi ? g_cfg.get< tag::fieldout_ >()[ m_meshid ]
1090 : 25068 : : g_cfg.get< tag::fieldout >();
1091 : :
1092 : 25068 : return m_it % ft.get< tag::iter >() == 0;
1093 : : }
1094 : :
1095 : : bool
1096 : 41081 : Discretization::fieldtime() const
1097 : : // *****************************************************************************
1098 : : // Decide if field output physics time interval is hit
1099 : : //! \return True if field output physics time interval is hit
1100 : : // *****************************************************************************
1101 : : {
1102 [ + + ]: 41081 : if (g_cfg.get< tag::benchmark >()) return false;
1103 : :
1104 : 22481 : auto eps = std::numeric_limits< tk::real >::epsilon();
1105 : 22481 : bool multi = g_cfg.get< tag::input >().size() > 1;
1106 [ - + ]: 22481 : const auto& ftab = multi ? g_cfg.get< tag::fieldout_ >()[ m_meshid ]
1107 : 22481 : : g_cfg.get< tag::fieldout >();
1108 : 22481 : auto ft = ftab.get< tag::time >();
1109 : :
1110 [ + + ]: 22481 : if (ft < eps) return false;
1111 : :
1112 : 22449 : return std::floor(m_t/ft) - m_physFieldFloor > eps;
1113 : : }
1114 : :
1115 : : bool
1116 : 41057 : Discretization::fieldrange() const
1117 : : // *****************************************************************************
1118 : : // Decide if physics time falls into a field output time range
1119 : : //! \return True if physics time falls into a field output time range
1120 : : // *****************************************************************************
1121 : : {
1122 [ + + ]: 41057 : if (g_cfg.get< tag::benchmark >()) return false;
1123 : :
1124 : 22457 : const auto eps = std::numeric_limits< tk::real >::epsilon();
1125 : :
1126 : 22457 : bool output = false;
1127 : :
1128 : 22457 : bool multi = g_cfg.get< tag::input >().size() > 1;
1129 [ - + ]: 22457 : const auto& ftab = multi ? g_cfg.get< tag::fieldout_ >()[ m_meshid ]
1130 : 22457 : : g_cfg.get< tag::fieldout >();
1131 : 22457 : const auto& rf = ftab.get< tag::range >();
1132 [ + + ]: 22957 : for (std::size_t i=0; i<rf.size(); ++i) {
1133 [ + + ][ + + ]: 500 : if (m_t > rf[i][0] and m_t < rf[i][1]) {
[ + + ]
1134 : 40 : output |= std::floor(m_t/rf[i][2]) - m_rangeFieldFloor[i] > eps;
1135 : : }
1136 : : }
1137 : :
1138 : 22457 : return output;
1139 : : }
1140 : :
1141 : : bool
1142 : 43668 : Discretization::histiter() const
1143 : : // *****************************************************************************
1144 : : // Decide if history output iteration count interval is hit
1145 : : //! \return True if history output iteration count interval is hit
1146 : : // *****************************************************************************
1147 : : {
1148 [ + + ]: 43668 : if (g_cfg.get< tag::benchmark >()) return false;
1149 : :
1150 : 25068 : auto hist = g_cfg.get< tag::histout, tag::iter >();
1151 : 25068 : const auto& hist_points = g_cfg.get< tag::histout, tag::points >();
1152 : :
1153 [ + + ][ + - ]: 25068 : return m_it % hist == 0 and not hist_points.empty();
1154 : : }
1155 : :
1156 : : bool
1157 : 42882 : Discretization::histtime() const
1158 : : // *****************************************************************************
1159 : : // Decide if history output physics time interval is hit
1160 : : //! \return True if history output physics time interval is hit
1161 : : // *****************************************************************************
1162 : : {
1163 [ + + ]: 42882 : if (g_cfg.get< tag::benchmark >()) return false;
1164 : :
1165 : 24282 : auto eps = std::numeric_limits< tk::real >::epsilon();
1166 : 24282 : auto ht = g_cfg.get< tag::histout, tag::time >();
1167 : :
1168 [ - + ]: 24282 : if (ht < eps) return false;
1169 : :
1170 : 24282 : return std::floor(m_t/ht) - m_physHistFloor > eps;
1171 : : }
1172 : :
1173 : : bool
1174 : 42624 : Discretization::histrange() const
1175 : : // *****************************************************************************
1176 : : // Decide if physics time falls into a history output time range
1177 : : //! \return True if physics time falls into a history output time range
1178 : : // *****************************************************************************
1179 : : {
1180 [ + + ]: 42624 : if (g_cfg.get< tag::benchmark >()) return false;
1181 : :
1182 : 24024 : auto eps = std::numeric_limits< tk::real >::epsilon();
1183 : :
1184 : 24024 : bool output = false;
1185 : :
1186 : 24024 : const auto& rh = g_cfg.get< tag::histout, tag::range >();
1187 [ + + ]: 26184 : for (std::size_t i=0; i<rh.size(); ++i) {
1188 [ + + ][ + + ]: 2160 : if (m_t > rh[i][0] and m_t < rh[i][1]) {
[ + + ]
1189 : 590 : output |= std::floor(m_t/rh[i][2]) - m_rangeHistFloor[i] > eps;
1190 : : }
1191 : : }
1192 : :
1193 : 24024 : return output;
1194 : : }
1195 : :
1196 : : bool
1197 : 43668 : Discretization::integiter() const
1198 : : // *****************************************************************************
1199 : : // Decide if integral output iteration count interval is hit
1200 : : //! \return True if integral output iteration count interval is hit
1201 : : // *****************************************************************************
1202 : : {
1203 [ + + ]: 43668 : if (g_cfg.get< tag::benchmark >()) return false;
1204 : :
1205 : 25068 : auto integ = g_cfg.get< tag::integout, tag::iter >();
1206 : 25068 : const auto& sidesets_integral = g_cfg.get< tag::integout, tag::sidesets >();
1207 : :
1208 [ + + ][ + - ]: 25068 : return m_it % integ == 0 and not sidesets_integral.empty();
1209 : : }
1210 : :
1211 : : bool
1212 : 43342 : Discretization::integtime() const
1213 : : // *****************************************************************************
1214 : : // Decide if integral output physics time interval is hit
1215 : : //! \return True if integral output physics time interval is hit
1216 : : // *****************************************************************************
1217 : : {
1218 [ + + ]: 43342 : if (g_cfg.get< tag::benchmark >()) return false;
1219 : :
1220 : 24742 : auto eps = std::numeric_limits< tk::real >::epsilon();
1221 : 24742 : auto it = g_cfg.get< tag::integout, tag::time >();
1222 : :
1223 [ - + ]: 24742 : if (it < eps) return false;
1224 : :
1225 : 24742 : return std::floor(m_t/it) - m_physIntegFloor > eps;
1226 : : }
1227 : :
1228 : : bool
1229 : 43315 : Discretization::integrange() const
1230 : : // *****************************************************************************
1231 : : // Decide if physics time falls into a integral output time range
1232 : : //! \return True if physics time falls into a integral output time range
1233 : : // *****************************************************************************
1234 : : {
1235 [ + + ]: 43315 : if (g_cfg.get< tag::benchmark >()) return false;
1236 : :
1237 : 24715 : auto eps = std::numeric_limits< tk::real >::epsilon();
1238 : :
1239 : 24715 : bool output = false;
1240 : :
1241 : 24715 : const auto& ri = g_cfg.get< tag::integout, tag::range >();
1242 [ + + ]: 24775 : for (std::size_t i=0; i<ri.size(); ++i) {
1243 [ + + ][ + + ]: 60 : if (m_t > ri[i][0] and m_t < ri[i][1]) {
[ + + ]
1244 : 15 : output |= std::floor(m_t/ri[i][2]) - m_rangeIntegFloor[i] > eps;
1245 : : }
1246 : : }
1247 : :
1248 : 24715 : return output;
1249 : : }
1250 : :
1251 : : bool
1252 : 45489 : Discretization::finished() const
1253 : : // *****************************************************************************
1254 : : // Decide if this is the last time step
1255 : : //! \return True if this is the last time step
1256 : : // *****************************************************************************
1257 : : {
1258 : 45489 : auto eps = std::numeric_limits< tk::real >::epsilon();
1259 : 45489 : auto nstep = g_cfg.get< tag::nstep >();
1260 : 45489 : auto term = g_cfg.get< tag::term >();
1261 : 45489 : auto residual = g_cfg.get< tag::residual >();
1262 : :
1263 [ + + ][ + + ]: 87979 : return std::abs(m_t-term) < eps or m_it >= nstep or
1264 [ + + ][ - + ]: 87979 : (m_res > 0.0 and m_res < residual);
1265 : : }
1266 : :
1267 : : void
1268 : 1060 : Discretization::residual( tk::real r )
1269 : : // *****************************************************************************
1270 : : // Update residual (during convergence to steady state)
1271 : : //! \param[in] r Current residual
1272 : : // *****************************************************************************
1273 : : {
1274 : 1060 : auto ttyi = g_cfg.get< tag::ttyi >();
1275 : :
1276 [ + - ]: 1060 : if (m_it % ttyi == 0) {
1277 : 1060 : m_res0 = m_res1;
1278 : 1060 : m_res1 = r;
1279 : : }
1280 : :
1281 : 1060 : m_res = r;
1282 : 1060 : }
1283 : :
1284 : : bool
1285 : 48 : Discretization::deastart()
1286 : : // *****************************************************************************
1287 : : // Decide whether to start the deactivation procedure in this time step
1288 : : //! \return True to start the deactivation prcedure in this time step
1289 : : // *****************************************************************************
1290 : : {
1291 [ + + ][ + - ]: 48 : if (not m_deastarted and m_t > g_cfg.get<tag::deatime>() + m_dt) {
[ + + ]
1292 : 19 : m_deastarted = 1;
1293 : 19 : return true;
1294 : : }
1295 : :
1296 : 29 : return false;
1297 : : }
1298 : :
1299 : : bool
1300 : 6086 : Discretization::deactivate() const
1301 : : // *****************************************************************************
1302 : : // Decide whether to run deactivation procedure in this time step
1303 : : //! \return True to run deactivation prcedure in this time step
1304 : : // *****************************************************************************
1305 : : {
1306 : 6086 : auto dea = g_cfg.get< tag::deactivate >();
1307 : 6086 : auto deafreq = g_cfg.get< tag::deafreq >();
1308 : :
1309 [ + + ][ + - ]: 6086 : if (dea and !m_nodeCommMap.empty() and m_it % deafreq == 0)
[ + + ][ + + ]
1310 : 200 : return true;
1311 : : else
1312 : 5886 : return false;
1313 : : }
1314 : :
1315 : : void
1316 : 10 : Discretization::deactivated( int n )
1317 : : // *****************************************************************************
1318 : : // Receive deactivation report
1319 : : //! \param[in] n Sum of deactivated chares
1320 : : // *****************************************************************************
1321 : : {
1322 : 10 : m_dea = n;
1323 : 10 : }
1324 : :
1325 : : bool
1326 : 17212 : Discretization::lb() const
1327 : : // *****************************************************************************
1328 : : // Decide whether to do load balancing this time step
1329 : : //! \return True to do load balancing in this time step
1330 : : // *****************************************************************************
1331 : : {
1332 : 17212 : auto lbfreq = g_cfg.get< tag::lbfreq >();
1333 : 17212 : auto lbtime = g_cfg.get< tag::lbtime >();
1334 : :
1335 [ + - ][ + + ]: 17212 : if ((m_t > lbtime and m_it % lbfreq == 0) or m_it == 2)
[ + + ]
1336 : 12064 : return true;
1337 : : else
1338 : 5148 : return false;
1339 : : }
1340 : :
1341 : : void
1342 : 701 : Discretization::pit( std::size_t it )
1343 : : // *****************************************************************************
1344 : : // Update number of pressure linear solve iterations taken
1345 : : //! \param[in] it Number of pressure linear solve iterations taken
1346 : : // *****************************************************************************
1347 : : {
1348 : 701 : m_pit = it;
1349 : 701 : }
1350 : :
1351 : : void
1352 : 20 : Discretization::mit( std::size_t it )
1353 : : // *****************************************************************************
1354 : : // Update number of momentum/transport linear solve iterations taken
1355 : : //! \param[in] it Number of momentum/transport linear solve iterations taken
1356 : : // *****************************************************************************
1357 : : {
1358 : 20 : m_mit = it;
1359 : 20 : }
1360 : :
1361 : : void
1362 : 255 : Discretization::npoin( std::size_t n )
1363 : : // *****************************************************************************
1364 : : // Set number of mesh points (across all meshes)
1365 : : //! \param[in] n Number of mesh points
1366 : : // *****************************************************************************
1367 : : {
1368 : 255 : m_npoin = n;
1369 : 255 : }
1370 : :
1371 : : void
1372 : 18577 : Discretization::status()
1373 : : // *****************************************************************************
1374 : : // Output one-liner status report
1375 : : // *****************************************************************************
1376 : : {
1377 : 18577 : const auto ttyi = g_cfg.get< tag::ttyi >();
1378 : :
1379 [ + + ][ + - ]: 18577 : if (thisIndex == 0 and m_meshid == 0 and m_it % ttyi == 0) {
[ + + ]
1380 : :
1381 : : // estimate grind time (taken between this and the previous time step)
1382 : : using std::chrono::duration_cast;
1383 : : using us = std::chrono::microseconds;
1384 : : using clock = std::chrono::high_resolution_clock;
1385 [ + - ][ + - ]: 3060 : auto grind_time = duration_cast< us >(clock::now() - m_prevstatus).count()
1386 : 3060 : / static_cast< long >( ttyi );
1387 : 3060 : auto grind_perf = static_cast<tk::real>(grind_time)
1388 : 3060 : / static_cast<tk::real>(m_npoin);
1389 : 3060 : m_prevstatus = clock::now();
1390 : :
1391 : 3060 : const auto term = g_cfg.get< tag::term >();
1392 : 3060 : const auto t0 = g_cfg.get< tag::t0 >();
1393 : 3060 : const auto nstep = g_cfg.get< tag::nstep >();
1394 : 3060 : const auto diag = g_cfg.get< tag::diag_iter >();
1395 : 3060 : const auto rsfreq = g_cfg.get< tag::rsfreq >();
1396 : 3060 : const auto benchmark = g_cfg.get< tag::benchmark >();
1397 : 3060 : const auto residual = g_cfg.get< tag::residual >();
1398 [ + - ]: 3060 : const auto solver = g_cfg.get< tag::solver >();
1399 [ + + ]: 3060 : const auto pre = solver == "chocg" ? 1 : 0;
1400 : 3060 : const auto theta = g_cfg.get< tag::theta >();
1401 : 3060 : const auto eps = std::numeric_limits< tk::real >::epsilon();
1402 [ + + ][ + + ]: 3060 : const auto mom = solver == "chocg" and theta > eps ? 1 : 0;
1403 : :
1404 : : // estimate time elapsed and time for accomplishment
1405 [ + - ][ + - ]: 3060 : tk::Timer::Watch ete, eta;
1406 [ + - ]: 3060 : m_timer.eta( term-t0, m_t-t0, nstep, m_it, m_res0, m_res1, residual,
1407 : : ete, eta );
1408 : :
1409 : 3060 : tk::Print print;
1410 : :
1411 : : // Output one-liner
1412 : 3060 : print << std::setfill(' ') << std::setw(8) << m_it << " "
1413 : 0 : << std::scientific << std::setprecision(6)
1414 : 0 : << std::setw(12) << m_t << " "
1415 : 3060 : << m_dt << " "
1416 : 0 : << std::setfill('0')
1417 : 0 : << std::setw(3) << ete.hrs.count() << ":"
1418 : 0 : << std::setw(2) << ete.min.count() << ":"
1419 : 0 : << std::setw(2) << ete.sec.count() << " "
1420 : 0 : << std::setw(3) << eta.hrs.count() << ":"
1421 : 0 : << std::setw(2) << eta.min.count() << ":"
1422 : 0 : << std::setw(2) << eta.sec.count() << " "
1423 : 0 : << std::scientific << std::setprecision(6) << std::setfill(' ')
1424 : 0 : << std::setw(9) << grind_time << " "
1425 [ + - ][ + - ]: 3060 : << std::setw(9) << grind_perf << " ";
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
1426 : :
1427 : : // Augment one-liner status with output indicators
1428 [ + - ][ + + ]: 3060 : if (fielditer() or fieldtime() or fieldrange()) print << 'f';
[ + - ][ + + ]
[ + - ][ + + ]
[ + + ][ + - ]
1429 [ + + ][ + - ]: 3060 : if (not (m_it % diag)) print << 'd';
1430 [ + - ][ + + ]: 3060 : if (histiter() or histtime() or histrange()) print << 't';
[ + - ][ + + ]
[ + - ][ + + ]
[ + + ][ + - ]
1431 [ + - ][ + + ]: 3060 : if (integiter() or integtime() or integrange()) print << 'i';
[ + - ][ + + ]
[ + - ][ + + ]
[ + + ][ + - ]
1432 [ - + ][ - - ]: 3060 : if (m_refined) print << 'h';
1433 [ + - ][ + + ]: 3060 : if (lb() and not finished()) print << 'l';
[ + - ][ + + ]
[ + + ][ + - ]
1434 [ + + ][ + - ]: 3060 : if (not benchmark && (not (m_it % rsfreq) || finished())) print << 'c';
[ + - ][ + + ]
[ + + ][ + - ]
1435 [ + + ][ + - ]: 3060 : if (m_deastarted and deactivate()) {
[ + + ][ + + ]
1436 [ + - ][ + - ]: 10 : print << "\te:" << m_dea << '/' << m_nchare;
[ + - ][ + - ]
1437 : : }
1438 [ + + ][ + - ]: 3060 : if (pre) print << "\tp:" << m_pit;
[ + - ]
1439 [ + + ][ + - ]: 3060 : if (mom) print << "\tm:" << m_mit;
[ + - ]
1440 : :
1441 [ + - ]: 3060 : print << '\n';
1442 : 3060 : }
1443 : 18577 : }
1444 : :
1445 : : #include "NoWarning/discretization.def.h"
|
