Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/Discretization.hpp
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 : : The Discretization class contains data and functionality common to all
11 : : discretization schemes.
12 : : */
13 : : // *****************************************************************************
14 : : #pragma once
15 : :
16 : : #include "Types.hpp"
17 : : #include "Timer.hpp"
18 : : #include "Fields.hpp"
19 : : #include "PUPUtil.hpp"
20 : : #include "UnsMesh.hpp"
21 : : #include "History.hpp"
22 : :
23 : : #include "NoWarning/discretization.decl.h"
24 : : #include "NoWarning/refiner.decl.h"
25 : :
26 : : namespace inciter {
27 : :
28 : : //! \brief Discretization Charm++ chare array holding common functinoality to
29 : : //! all discretization schemes
30 : : class Discretization : public CBase_Discretization {
31 : :
32 : : public:
33 : : #if defined(__clang__)
34 : : #pragma clang diagnostic push
35 : : #pragma clang diagnostic ignored "-Wunused-parameter"
36 : : #pragma clang diagnostic ignored "-Wdeprecated-declarations"
37 : : #elif defined(STRICT_GNUC)
38 : : #pragma GCC diagnostic push
39 : : #pragma GCC diagnostic ignored "-Wunused-parameter"
40 : : #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
41 : : #elif defined(__INTEL_COMPILER)
42 : : #pragma warning( push )
43 : : #pragma warning( disable: 1478 )
44 : : #endif
45 : : // Include Charm++ SDAG code. See http://charm.cs.illinois.edu/manuals/html/
46 : : // charm++/manual.html, Sec. "Structured Control Flow: Structured Dagger".
47 : : Discretization_SDAG_CODE
48 : : #if defined(__clang__)
49 : : #pragma clang diagnostic pop
50 : : #elif defined(STRICT_GNUC)
51 : : #pragma GCC diagnostic pop
52 : : #elif defined(__INTEL_COMPILER)
53 : : #pragma warning( pop )
54 : : #endif
55 : :
56 : : //! Constructor
57 : : explicit
58 : : Discretization(
59 : : std::size_t meshid,
60 : : const std::vector< CProxy_Discretization >& disc,
61 : : const CProxy_Transporter& transporter,
62 : : const tk::CProxy_MeshWriter& meshwriter,
63 : : const tk::UnsMesh::CoordMap& coordmap,
64 : : const tk::UnsMesh::Chunk& el,
65 : : const std::map< int, std::unordered_set< std::size_t > >& nodeCommMap,
66 : : int nc );
67 : :
68 : : #if defined(__clang__)
69 : : #pragma clang diagnostic push
70 : : #pragma clang diagnostic ignored "-Wundefined-func-template"
71 : : #endif
72 : : //! Migrate constructor
73 : : // cppcheck-suppress uninitMemberVar
74 : 4422 : explicit Discretization( CkMigrateMessage* m ) : CBase_Discretization( m )
75 : 4422 : {}
76 : : #if defined(__clang__)
77 : : #pragma clang diagnostic pop
78 : : #endif
79 : :
80 : : //! Configure Charm++ reduction types
81 : : static void registerReducers();
82 : :
83 : : //! Resize mesh data structures after mesh refinement
84 : : void resizePostAMR(
85 : : const tk::UnsMesh::Chunk& chunk,
86 : : const tk::UnsMesh::Coords& coord,
87 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
88 : : nodeCommMap,
89 : : const std::set< std::size_t >& removedNodes );
90 : :
91 : : //! Get ready for (re-)computing/communicating nodal volumes
92 : : void startvol();
93 : :
94 : : //! Sum mesh volumes to nodes, start communicating them on chare-boundaries
95 : : void vol();
96 : :
97 : : //! Set Refiner Charm++ proxy
98 : : void setRefiner( const CProxy_Refiner& ref );
99 : :
100 : : //! Collect nodal volumes across chare boundaries
101 : : void comvol( int c,
102 : : const std::vector< std::size_t >& gid,
103 : : const std::vector< tk::real >& nodevol );
104 : :
105 : : //! Sum mesh volumes and contribute own mesh volume to total volume
106 : : void totalvol();
107 : :
108 : : //! Compute mesh cell statistics
109 : : void stat( tk::real mesh_volume );
110 : :
111 : : //! Decide whether to start the deactivation procedure in this time step
112 : : bool deastart();
113 : :
114 : : //! Decide whether to run deactivation procedure in this time step
115 : : bool deactivate() const;
116 : :
117 : : //! Receive deactivation report
118 : : void deactivated( int n );
119 : :
120 : : //! Decide whether to do load balancing this time step
121 : : bool lb() const;
122 : :
123 : : //! Compute total box IC volume
124 : : void boxvol();
125 : :
126 : : /** @name Accessors */
127 : : ///@{
128 : : //! Coordinates accessor as const-ref
129 : : const tk::UnsMesh::Coords& Coord() const { return m_coord; }
130 : : //! Coordinates accessor as reference
131 : 306556 : tk::UnsMesh::Coords& Coord() { return m_coord; }
132 : :
133 : : //! Global ids accessors as const-ref
134 : 14467 : const std::vector< std::size_t >& Gid() const { return m_gid; }
135 : :
136 : : //! Local ids accessors as const-ref
137 : 744964 : const std::unordered_map< std::size_t, std::size_t >& Lid() const
138 : 744964 : { return m_lid; }
139 : :
140 : : //! Tetrahedron element connectivity (with local ids) accessors as const-ref
141 : 36484 : const std::vector< std::size_t >& Inpoel() const { return m_inpoel; }
142 : :
143 : : //! Mesh chunk accessor as const-ref
144 : : const tk::UnsMesh::Chunk& Chunk() const { return m_el; }
145 : :
146 : : //! Total mesh volume accessor
147 : : tk::real meshvol() const { return m_meshvol; }
148 : :
149 : : //! Nodal mesh volume accessors const-ref
150 : 103250 : const std::vector< tk::real >& V() const { return m_v; }
151 : : //! Nodal mesh volumes at current time step accessors as const-ref
152 : 171260 : const std::vector< tk::real >& Vol() const { return m_vol; }
153 : : //! Access chare-volume contributions per chare per chare-boundary node
154 : : const std::unordered_map< int,
155 : 259 : std::unordered_map< std::size_t, tk::real > >& Cvolc() const {
156 : 259 : return m_cvolc;
157 : : }
158 : :
159 : : //! Access number of chares
160 : 190 : int nchare() const { return m_nchare; }
161 : :
162 : : //! Set 'initial' flag
163 : : //! \param[in] i Value to put in 'initial'
164 : 2544 : void Initial( std::size_t i ) { m_initial = i; }
165 : : //! Query 'initial' flag
166 : : //! \return True during setup, false durign time stepping
167 : 16006 : bool Initial() const { return m_initial; }
168 : :
169 : : //! History points data accessor as const-ref
170 : 2476 : const std::vector< HistData >& Hist() const { return m_histdata; }
171 : :
172 : : //! Box volume accessor
173 : 2576 : tk::real& Boxvol() { return m_boxvol; }
174 : :
175 : : //! Mesh ID accessor
176 : 139420 : std::size_t MeshId() const { return m_meshid; }
177 : :
178 : : //! Time step size accessor
179 : 1906378 : tk::real Dt() const { return m_dt; }
180 : : //! Time step size at previous time step accessor
181 : : tk::real Dtn() const { return m_dtn; }
182 : : //! Physical time accessor
183 : 1583796 : tk::real T() const { return m_t; }
184 : : //! Iteration count accessor
185 : 159569 : uint64_t It() const { return m_it; }
186 : :
187 : : //! Non-const-ref refinement iteration count accessor
188 : 0 : uint64_t& Itr() { return m_itr; }
189 : : //! Non-const-ref field-output iteration count accessor
190 : 0 : uint64_t& Itf() { return m_itf; }
191 : :
192 : : //! Non-const-ref number of restarts accessor
193 : : int& Nrestart() { return m_nrestart; }
194 : :
195 : : //! Timer accessor as const-ref
196 : : const tk::Timer& Timer() const { return m_timer; }
197 : : //! Timer accessor as non-const-ref
198 : 2544 : tk::Timer& Timer() { return m_timer; }
199 : :
200 : : //! Accessor to flag indicating if the mesh was refined as a value
201 : : int refined() const { return m_refined; }
202 : : //! Accessor to flag indicating if the mesh was refined as non-const-ref
203 : 40608 : int& refined() { return m_refined; }
204 : :
205 : : //! Transporter proxy accessor as const-ref
206 : : const CProxy_Transporter& Tr() const { return m_transporter; }
207 : : //! Transporter proxy accessor as non-const-ref
208 : 46936 : CProxy_Transporter& Tr() { return m_transporter; }
209 : :
210 : : //! Access bound Refiner class pointer
211 : 0 : Refiner* Ref() const {
212 [ - - ][ - - ]: 0 : Assert( m_refiner[ thisIndex ].ckLocal() != nullptr,
[ - - ][ - - ]
[ - - ][ - - ]
213 : : "Refiner ckLocal() null" );
214 [ - - ][ - - ]: 0 : return m_refiner[ thisIndex ].ckLocal();
215 : : }
216 : :
217 : : //! Node communication map accessor as const-ref
218 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
219 : 1986692 : NodeCommMap() const { return m_nodeCommMap; }
220 : :
221 : : //! IC boxnodes accessor
222 : 97480 : const std::vector< std::unordered_set< std::size_t > >& BoxNodes() const
223 : 97480 : { return m_boxnodes; }
224 : :
225 : : //! Transfer flags accessor
226 : 0 : const std::vector< double >& TransferFlag() const { return m_transfer_flag; }
227 : : //@}
228 : :
229 : : //! Set time step size
230 : : void setdt( tk::real newdt );
231 : :
232 : : //! Prepare for next step
233 : : void next();
234 : :
235 : : //! Otput one-liner status report
236 : : void status();
237 : :
238 : : //! Construct history output filename
239 : : std::string histfilename( const std::string& id,
240 : : std::streamsize precision );
241 : :
242 : : //! Output headers for time history files (one for each point)
243 : : void histheader( std::vector< std::string >&& names );
244 : :
245 : : //! Output time history for a time step
246 : : void history( std::vector< std::vector< tk::real > >&& data );
247 : :
248 : : //! Our mesh has been registered with mesh-to-mesh transfer (if coupled)
249 : : void transfer_initialized();
250 : :
251 : : //! Initiate solution transfer from background to overset mesh (if coupled)
252 : : void transfer( tk::Fields& u, CkCallback c, bool trflag );
253 : :
254 : : //! Initiate solution transfer from overset to background mesh
255 : : void transfer_from();
256 : :
257 : : //! Prepare integrid-boundary data structures (if coupled)
258 : : void intergrid( const std::map< int, std::vector< std::size_t > >& bnode );
259 : :
260 : : //! Communicate holes to background mesh
261 : : void hole( const std::map< int, std::vector< std::size_t > >& bface,
262 : : const std::vector< std::size_t >& triinpoel,
263 : : CkCallback c );
264 : :
265 : : //! Receive hole data
266 : : void aggregateHoles( CkReductionMsg* msg );
267 : :
268 : : //! Hole communication complete
269 : : void holeComplete();
270 : :
271 : : //! Find mesh nodes within hols
272 : : void holefind();
273 : :
274 : : //! Output mesh and fields data (solution dump) to file(s)
275 : : void write( const std::vector< std::size_t >& inpoel,
276 : : const tk::UnsMesh::Coords& coord,
277 : : const std::map< int, std::vector< std::size_t > >& bface,
278 : : const std::map< int, std::vector< std::size_t > >& bnode,
279 : : const std::vector< std::size_t >& triinpoel,
280 : : const std::vector< std::string>& elemfieldnames,
281 : : const std::vector< std::string>& nodefieldnames,
282 : : const std::vector< std::string>& elemsurfnames,
283 : : const std::vector< std::string>& nodesurfnames,
284 : : const std::vector< std::vector< tk::real > >& elemfields,
285 : : const std::vector< std::vector< tk::real > >& nodefields,
286 : : const std::vector< std::vector< tk::real > >& elemsurfs,
287 : : const std::vector< std::vector< tk::real > >& nodesurfs,
288 : : CkCallback c );
289 : :
290 : : //! Zero grind-timer
291 : : void grindZero();
292 : :
293 : : //! Detect if just returned from a checkpoint and if so, zero timers
294 : : bool restarted( int nrestart );
295 : :
296 : : //! Remap mesh data due to new local ids
297 : : void remap( const std::unordered_map< std::size_t, std::size_t >& map );
298 : :
299 : : //! Decide if field output iteration count interval is hit
300 : : bool fielditer() const;
301 : : //! Decide if field output physics time interval is hit
302 : : bool fieldtime() const;
303 : : //! Decide if physics time falls into a field output time range
304 : : bool fieldrange() const;
305 : :
306 : : //! Decide if history output iteration count interval is hit
307 : : bool histiter() const;
308 : : //! Decide if history output physics time interval is hit
309 : : bool histtime() const;
310 : : //! Decide if physics time falls into a history output time range
311 : : bool histrange() const;
312 : :
313 : : //! Decide if integral output iteration count interval is hit
314 : : bool integiter() const;
315 : : //! Decide if integral output physics time interval is hit
316 : : bool integtime() const;
317 : : //! Decide if physics time falls into a integral output time range
318 : : bool integrange() const;
319 : :
320 : : //! Decide if this is the last time step
321 : : bool finished() const;
322 : :
323 : : //! Update residual (during convergence to steady state)
324 : : void residual( tk::real r );
325 : :
326 : : //! Update number of pressure linear solve iterations taken
327 : : void pit( std::size_t it );
328 : :
329 : : //! Update number of momentum/transport linear solve iterations taken
330 : : void mit( std::size_t it );
331 : :
332 : : //! Set number of mesh points (across all meshes)
333 : : void npoin( std::size_t n );
334 : :
335 : : /** @name Charm++ pack/unpack serializer member functions */
336 : : ///@{
337 : : //! \brief Pack/Unpack serialize member function
338 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
339 : 14060 : void pup( PUP::er &p ) override {
340 : 14060 : p | m_npoin;
341 : 14060 : p | m_meshid;
342 : 14060 : p | m_transfer_complete;
343 : 14060 : p | m_transfer_flag;
344 : 14060 : p | m_trflag;
345 : 14060 : p | m_transfer_hole;
346 : 14060 : p | m_holcont;
347 : 14060 : p | m_nchare;
348 : 14060 : p | m_it;
349 : 14060 : p | m_itr;
350 : 14060 : p | m_itf;
351 : 14060 : p | m_initial;
352 : 14060 : p | m_t;
353 : 14060 : p | m_lastDumpTime;
354 : 14060 : p | m_physFieldFloor;
355 : 14060 : p | m_physHistFloor;
356 : 14060 : p | m_physIntegFloor;
357 : 14060 : p | m_rangeFieldFloor;
358 : 14060 : p | m_rangeHistFloor;
359 : 14060 : p | m_rangeIntegFloor;
360 : 14060 : p | m_dt;
361 : 14060 : p | m_dtn;
362 : 14060 : p | m_nvol;
363 : 14060 : p | m_boxnodes;
364 : 14060 : p | m_disc;
365 : 14060 : p | m_transporter;
366 : 14060 : p | m_meshwriter;
367 : 14060 : p | m_refiner;
368 : 14060 : p | m_el;
369 [ + + ]: 14060 : if (p.isUnpacking()) {
370 : 4422 : m_inpoel = std::get< 0 >( m_el );
371 : 4422 : m_gid = std::get< 1 >( m_el );
372 : 4422 : m_lid = std::get< 2 >( m_el );
373 : : }
374 : 14060 : p | m_coord;
375 : 14060 : p | m_nodeCommMap;
376 : 14060 : p | m_meshvol;
377 : 14060 : p | m_v;
378 : 14060 : p | m_vol;
379 : 14060 : p | m_volc;
380 : 14060 : p | m_cvolc;
381 : 14060 : p | m_boxvol;
382 : 14060 : p | m_timer;
383 : 14060 : p | m_refined;
384 : 14060 : p( reinterpret_cast<char*>(&m_prevstatus), sizeof(Clock::time_point) );
385 : 14060 : p | m_nrestart;
386 : 14060 : p | m_histdata;
387 : 14060 : p | m_res;
388 : 14060 : p | m_res0;
389 : 14060 : p | m_res1;
390 : 14060 : p | m_dea;
391 : 14060 : p | m_deastarted;
392 : 14060 : p | m_pit;
393 : 14060 : p | m_mit;
394 : 14060 : }
395 : : //! \brief Pack/Unpack serialize operator|
396 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
397 : : //! \param[in,out] i Discretization object reference
398 : : friend void operator|( PUP::er& p, Discretization& i ) { i.pup(p); }
399 : : //@}
400 : :
401 : : private:
402 : : // Shorthand for clock, setting an internal clock type
403 : : using Clock = std::chrono::high_resolution_clock;
404 : :
405 : : //! Total number of mesh points (across all meshes)
406 : : std::size_t m_npoin;
407 : : //! Mesh ID
408 : : std::size_t m_meshid;
409 : : //! Function to call after mesh-to-mesh solution transfer is complete
410 : : CkCallback m_transfer_complete;
411 : : //! Pointer to solution during mesh-to-mesh solution transfer
412 : : tk::Fields* m_transfer_sol;
413 : : //! Transfer flags at nodes (if coupled)
414 : : std::vector< double > m_transfer_flag;
415 : : //! Transfer flags if true
416 : : bool m_trflag;
417 : : //! Holes tessellation
418 : : //! \details: Key: hole id, value: list of triangle node coordinates
419 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_transfer_hole;
420 : : //! Callback to continue with after hole communication
421 : : CkCallback m_holcont;
422 : : //! Total number of Discretization chares
423 : : int m_nchare;
424 : : //! Iteration count
425 : : uint64_t m_it;
426 : : //! Iteration count with mesh refinement
427 : : //! \details Used as the restart sequence number {RS} in saving output in
428 : : //! an ExodusII sequence
429 : : //! \see https://www.paraview.org/Wiki/Restarted_Simulation_Readers
430 : : uint64_t m_itr;
431 : : //! Field output iteration count without mesh refinement
432 : : //! \details Counts the number of field outputs to file during two
433 : : //! time steps with mesh efinement
434 : : uint64_t m_itf;
435 : : //! Flag that is nonzero during setup and zero during time stepping
436 : : std::size_t m_initial;
437 : : //! Physical time
438 : : tk::real m_t;
439 : : //! Physics time at last field output
440 : : tk::real m_lastDumpTime;
441 : : //! Recent floor of physics time divided by field output interval time
442 : : tk::real m_physFieldFloor;
443 : : //! Recent floor of physics time divided by history output interval time
444 : : tk::real m_physHistFloor;
445 : : //! Recent floor of physics time divided by integral output interval time
446 : : tk::real m_physIntegFloor;
447 : : //! Recent floors of physics time divided by field output time for ranges
448 : : std::vector< tk::real > m_rangeFieldFloor;
449 : : //! Recent floors of physics time divided by history output time for ranges
450 : : std::vector< tk::real > m_rangeHistFloor;
451 : : //! Recent floors of physics time divided by integral output time for ranges
452 : : std::vector< tk::real > m_rangeIntegFloor;
453 : : //! Physical time step size
454 : : tk::real m_dt;
455 : : //! Physical time step size at the previous time step
456 : : tk::real m_dtn;
457 : : //! \brief Number of chares from which we received nodal volume
458 : : //! contributions on chare boundaries
459 : : std::size_t m_nvol;
460 : : //! List of nodes at which box user ICs are set for each IC box
461 : : std::vector< std::unordered_set< std::size_t > > m_boxnodes;
462 : : //! Discretization proxy for all meshes
463 : : std::vector< CProxy_Discretization > m_disc;
464 : : //! Transporter proxy
465 : : CProxy_Transporter m_transporter;
466 : : //! Mesh writer proxy
467 : : tk::CProxy_MeshWriter m_meshwriter;
468 : : //! Mesh refiner proxy
469 : : CProxy_Refiner m_refiner;
470 : : //! \brief Elements of the mesh chunk we operate on
471 : : //! \details Initialized by the constructor. The first vector is the element
472 : : //! connectivity (local IDs), the second vector is the global node IDs of
473 : : //! owned elements, while the third one is a map of global->local node
474 : : //! IDs.
475 : : tk::UnsMesh::Chunk m_el;
476 : : //! Alias to element connectivity
477 : : std::vector< std::size_t >& m_inpoel = std::get<0>( m_el );
478 : : //! Alias to global node IDs of owned elements
479 : : std::vector< std::size_t >& m_gid = std::get<1>( m_el );
480 : : //! \brief Alias to local node ids associated to the global ones of owned
481 : : //! elements
482 : : std::unordered_map< std::size_t, std::size_t >& m_lid = std::get<2>( m_el );
483 : : //! Mesh point coordinates
484 : : tk::UnsMesh::Coords m_coord;
485 : : //! \brief Global mesh node IDs bordering the mesh chunk held by fellow
486 : : //! Discretization chares associated to their chare IDs
487 : : std::unordered_map< int, std::unordered_set< std::size_t > > m_nodeCommMap;
488 : : //! Total mesh volume
489 : : tk::real m_meshvol;
490 : : //! Nodal mesh volumes
491 : : //! \details This is the volume of the mesh associated to nodes of owned
492 : : //! elements (sum of surrounding cell volumes / 4) without contributions
493 : : //! from other chares on chare-boundaries
494 : : std::vector< tk::real > m_v;
495 : : //! Volume of nodes
496 : : //! \details This is the volume of the mesh associated to nodes of owned
497 : : //! elements (sum of surrounding cell volumes / 4) with contributions from
498 : : //! other chares on chare-boundaries
499 : : std::vector< tk::real > m_vol;
500 : : //! Receive buffer for volume of nodes (with global node id as key)
501 : : //! \details This is a communication buffer used to compute the volume of
502 : : //! the mesh associated to nodes of owned elements (sum of surrounding
503 : : //! cell volumes / 4) with contributions from other chares on
504 : : //! chare-boundaries.
505 : : std::unordered_map< std::size_t, tk::real > m_volc;
506 : : //! Chare-volume contributions per chare per chare-boundary node
507 : : //! \details Outer key: chare id, value: local node id and node volume
508 : : std::unordered_map< int, std::unordered_map< std::size_t, tk::real > >
509 : : m_cvolc;
510 : : //! Volume of user-defined box IC
511 : : tk::real m_boxvol;
512 : : //! Timer measuring a time step
513 : : tk::Timer m_timer;
514 : : //! 1 if mesh was refined in a time step, 0 if it was not
515 : : int m_refined;
516 : : //! Time point storing clock state at status()
517 : : Clock::time_point m_prevstatus;
518 : : //! Number of times restarted
519 : : int m_nrestart;
520 : : //! Data at history point locations
521 : : std::vector< HistData > m_histdata;
522 : : //! Current residual (during convergence to steady state)
523 : : tk::real m_res;
524 : : //! Residual at previous ETA calcuation (during convergence to steady state)
525 : : tk::real m_res0;
526 : : //! Residual at next ETA calcuation (during convergence to steady state)
527 : : tk::real m_res1;
528 : : //! Numberf of deactived chares
529 : : int m_dea;
530 : : //! Flag: 1 if deactivation procedure has started
531 : : int m_deastarted;
532 : : //! Numberf of recent pressure solve iterations taken
533 : : std::size_t m_pit;
534 : : //! Numberf of recent momentum/transport solve iterations taken
535 : : std::size_t m_mit;
536 : :
537 : : //! Set mesh coordinates based on coordinates map
538 : : tk::UnsMesh::Coords setCoord( const tk::UnsMesh::CoordMap& coordmap );
539 : : };
540 : :
541 : : } // inciter::
|
