1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
// *****************************************************************************
/*!
  \file      src/Inciter/Discretization.hpp
  \copyright 2012-2015 J. Bakosi,
             2016-2018 Los Alamos National Security, LLC.,
             2019-2021 Triad National Security, LLC.,
             2022-2025 J. Bakosi
             All rights reserved. See the LICENSE file for details.
  \details   Data and functionality common to all discretization schemes
     The Discretization class contains data and functionality common to all
     discretization schemes.
*/
// *****************************************************************************
#pragma once

#include "Types.hpp"
#include "Timer.hpp"
#include "Fields.hpp"
#include "PUPUtil.hpp"
#include "UnsMesh.hpp"<--- Include file: "UnsMesh.hpp" not found.
#include "History.hpp"

#include "NoWarning/discretization.decl.h"<--- Include file: "NoWarning/discretization.decl.h" not found.
#include "NoWarning/refiner.decl.h"<--- Include file: "NoWarning/refiner.decl.h" not found.

namespace inciter {

//! \brief Discretization Charm++ chare array holding common functinoality to
//!   all discretization schemes
class Discretization : public CBase_Discretization {

  public:
    #if defined(__clang__)
      #pragma clang diagnostic push
      #pragma clang diagnostic ignored "-Wunused-parameter"
      #pragma clang diagnostic ignored "-Wdeprecated-declarations"
    #elif defined(STRICT_GNUC)
      #pragma GCC diagnostic push
      #pragma GCC diagnostic ignored "-Wunused-parameter"
      #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
    #elif defined(__INTEL_COMPILER)
      #pragma warning( push )
      #pragma warning( disable: 1478 )
    #endif
    // Include Charm++ SDAG code. See http://charm.cs.illinois.edu/manuals/html/
    // charm++/manual.html, Sec. "Structured Control Flow: Structured Dagger".
    Discretization_SDAG_CODE
    #if defined(__clang__)
      #pragma clang diagnostic pop
    #elif defined(STRICT_GNUC)
      #pragma GCC diagnostic pop
    #elif defined(__INTEL_COMPILER)
      #pragma warning( pop )
    #endif

    //! Constructor
    explicit
      Discretization(
        std::size_t meshid,
        const std::vector< CProxy_Discretization >& disc,
        const CProxy_Transporter& transporter,
        const tk::CProxy_MeshWriter& meshwriter,
        const tk::UnsMesh::CoordMap& coordmap,
        const tk::UnsMesh::Chunk& el,
        const std::map< int, std::unordered_set< std::size_t > >& nodeCommMap,
        int nc );

    #if defined(__clang__)
      #pragma clang diagnostic push
      #pragma clang diagnostic ignored "-Wundefined-func-template"
    #endif
    //! Migrate constructor
    // cppcheck-suppress uninitMemberVar
    explicit Discretization( CkMigrateMessage* m ) : CBase_Discretization( m )
    {}
    #if defined(__clang__)
      #pragma clang diagnostic pop
    #endif

    //! Configure Charm++ reduction types
    static void registerReducers();

    //! Resize mesh data structures after mesh refinement
    void resizePostAMR(
      const tk::UnsMesh::Chunk& chunk,
      const tk::UnsMesh::Coords& coord,
      const std::unordered_map< int, std::unordered_set< std::size_t > >&
        nodeCommMap,
      const std::set< std::size_t >& removedNodes );

    //! Get ready for (re-)computing/communicating nodal volumes
    void startvol();

    //! Sum mesh volumes to nodes, start communicating them on chare-boundaries
    void vol();

    //! Set Refiner Charm++ proxy
    void setRefiner( const CProxy_Refiner& ref );

    //! Collect nodal volumes across chare boundaries
    void comvol( int c,
                 const std::vector< std::size_t >& gid,
                 const std::vector< tk::real >& nodevol );

    //! Sum mesh volumes and contribute own mesh volume to total volume
    void totalvol();

    //! Compute mesh cell statistics
    void stat( tk::real mesh_volume );

    //! Decide whether to start the deactivation procedure in this time step
    bool deastart();

    //! Decide whether to run deactivation procedure in this time step
    bool deactivate() const;

    //!  Receive deactivation report
    void deactivated( int n );

    //! Decide whether to do load balancing this time step
    bool lb() const;

    //! Compute total box IC volume
    void boxvol();

    /** @name Accessors */
    ///@{
    //! Coordinates accessor as const-ref
    const tk::UnsMesh::Coords& Coord() const { return m_coord; }
    //! Coordinates accessor as reference
    tk::UnsMesh::Coords& Coord() { return m_coord; }

    //! Global ids accessors as const-ref
    const std::vector< std::size_t >& Gid() const { return m_gid; }

    //! Local ids accessors as const-ref
     const std::unordered_map< std::size_t, std::size_t >& Lid() const
    { return m_lid; }

    //! Tetrahedron element connectivity (with local ids) accessors as const-ref
    const std::vector< std::size_t >& Inpoel() const { return m_inpoel; }

    //! Mesh chunk accessor as const-ref
    const tk::UnsMesh::Chunk& Chunk() const { return m_el; }<--- The function 'Chunk' is never used.

    //! Total mesh volume accessor
    tk::real meshvol() const { return m_meshvol; }<--- The function 'meshvol' is never used.

    //! Nodal mesh volume accessors const-ref
    const std::vector< tk::real >& V() const { return m_v; }
    //! Nodal mesh volumes at current time step accessors as const-ref
    const std::vector< tk::real >& Vol() const { return m_vol; }
    //! Access chare-volume contributions per chare per chare-boundary node
    const std::unordered_map< int,
            std::unordered_map< std::size_t, tk::real > >& Cvolc() const {
      return m_cvolc;
    }

    //! Access number of chares
    int nchare() const { return m_nchare; }

    //! Set 'initial' flag
    //! \param[in] i Value to put in 'initial'
    void Initial( std::size_t i ) { m_initial = i; }
    //! Query 'initial' flag
    //! \return True during setup, false durign time stepping
    bool Initial() const { return m_initial; }

    //! History points data accessor as const-ref
    const std::vector< HistData >& Hist() const { return m_histdata; }

    //! Box volume accessor
    tk::real& Boxvol() { return m_boxvol; }

    //! Mesh ID accessor
    std::size_t MeshId() const { return m_meshid; }

    //! Time step size accessor
    tk::real Dt() const { return m_dt; }
    //! Time step size at previous time step accessor
    tk::real Dtn() const { return m_dtn; }<--- The function 'Dtn' is never used.
    //! Physical time accessor
    tk::real T() const { return m_t; }
    //! Iteration count accessor
    uint64_t It() const { return m_it; }

    //! Non-const-ref refinement iteration count accessor
    uint64_t& Itr() { return m_itr; }
    //! Non-const-ref field-output iteration count accessor
    uint64_t& Itf() { return m_itf; }

    //! Non-const-ref number of restarts accessor
    int& Nrestart() { return m_nrestart; }<--- The function 'Nrestart' is never used.

    //! Timer accessor as const-ref
    const tk::Timer& Timer() const { return m_timer; }
    //! Timer accessor as non-const-ref
    tk::Timer& Timer() { return m_timer; }

    //! Accessor to flag indicating if the mesh was refined as a value
    int refined() const { return m_refined; }
    //! Accessor to flag indicating if the mesh was refined as non-const-ref
    int& refined() { return m_refined; }

    //! Transporter proxy accessor as const-ref
    const CProxy_Transporter& Tr() const { return m_transporter; }
    //! Transporter proxy accessor as non-const-ref
    CProxy_Transporter& Tr() { return m_transporter; }

    //! Access bound Refiner class pointer
    Refiner* Ref() const {
      Assert( m_refiner[ thisIndex ].ckLocal() != nullptr,
              "Refiner ckLocal() null" );
      return m_refiner[ thisIndex ].ckLocal();
    }

    //! Node communication map accessor as const-ref
    const std::unordered_map< int, std::unordered_set< std::size_t > >&
      NodeCommMap() const { return m_nodeCommMap; }

    //! IC boxnodes accessor
    const std::vector< std::unordered_set< std::size_t > >& BoxNodes() const
    { return m_boxnodes; }

    //! Transfer flags accessor
    const std::vector< double >& TransferFlag() const { return m_transfer_flag; }
    //@}

    //! Set time step size
    void setdt( tk::real newdt );

    //! Prepare for next step
    void next();

    //! Otput one-liner status report
    void status();

    //! Construct history output filename
    std::string histfilename( const std::string& id,
                              std::streamsize precision );

    //! Output headers for time history files (one for each point)
    void histheader( std::vector< std::string >&& names );

    //! Output time history for a time step
    void history( std::vector< std::vector< tk::real > >&& data );

    //! Our mesh has been registered with mesh-to-mesh transfer (if coupled)
    void transfer_initialized();

    //! Initiate solution transfer from background to overset mesh (if coupled)
    void transfer( tk::Fields& u, CkCallback c, bool trflag );

    //! Initiate solution transfer from overset to background mesh
    void transfer_from();

    //! Prepare integrid-boundary data structures (if coupled)
    void intergrid( const std::map< int, std::vector< std::size_t > >& bnode );

    //! Communicate holes to background mesh
    void hole( const std::map< int, std::vector< std::size_t > >& bface,
               const std::vector< std::size_t >& triinpoel,
               CkCallback c );

    //! Receive hole data
    void aggregateHoles( CkReductionMsg* msg );

    //! Hole communication complete
    void holeComplete();

    //! Find mesh nodes within hols
    void holefind();

    //! Output mesh and fields data (solution dump) to file(s)
    void write( const std::vector< std::size_t >& inpoel,
                const tk::UnsMesh::Coords& coord,
                const std::map< int, std::vector< std::size_t > >& bface,
                const std::map< int, std::vector< std::size_t > >& bnode,
                const std::vector< std::size_t >& triinpoel,
                const std::vector< std::string>& elemfieldnames,
                const std::vector< std::string>& nodefieldnames,
                const std::vector< std::string>& elemsurfnames,
                const std::vector< std::string>& nodesurfnames,
                const std::vector< std::vector< tk::real > >& elemfields,
                const std::vector< std::vector< tk::real > >& nodefields,
                const std::vector< std::vector< tk::real > >& elemsurfs,
                const std::vector< std::vector< tk::real > >& nodesurfs,
                CkCallback c );

    //! Zero grind-timer
    void grindZero();

    //! Detect if just returned from a checkpoint and if so, zero timers
    bool restarted( int nrestart );

    //! Remap mesh data due to new local ids
    void remap( const std::unordered_map< std::size_t, std::size_t >& map );

    //! Decide if field output iteration count interval is hit
    bool fielditer() const;
    //! Decide if field output physics time interval is hit
    bool fieldtime() const;
    //! Decide if physics time falls into a field output time range
    bool fieldrange() const;

    //! Decide if history output iteration count interval is hit
    bool histiter() const;
    //! Decide if history output physics time interval is hit
    bool histtime() const;
    //! Decide if physics time falls into a history output time range
    bool histrange() const;

    //! Decide if integral output iteration count interval is hit
    bool integiter() const;
    //! Decide if integral output physics time interval is hit
    bool integtime() const;
    //! Decide if physics time falls into a integral output time range
    bool integrange() const;

    //! Decide if this is the last time step
    bool finished() const;

    //! Update residual (during convergence to steady state)
    void residual( tk::real r );

    //! Update number of pressure linear solve iterations taken
    void pit( std::size_t it );

    //! Update number of momentum/transport linear solve iterations taken
    void mit( std::size_t it );

    //! Set number of mesh points (across all meshes)
    void npoin( std::size_t n );

    /** @name Charm++ pack/unpack serializer member functions */
    ///@{
    //! \brief Pack/Unpack serialize member function
    //! \param[in,out] p Charm++'s PUP::er serializer object reference
    void pup( PUP::er &p ) override {
      p | m_npoin;
      p | m_meshid;
      p | m_transfer_complete;
      p | m_transfer_flag;
      p | m_trflag;
      p | m_transfer_hole;
      p | m_holcont;
      p | m_nchare;
      p | m_it;
      p | m_itr;
      p | m_itf;
      p | m_initial;
      p | m_t;
      p | m_lastDumpTime;
      p | m_physFieldFloor;
      p | m_physHistFloor;
      p | m_physIntegFloor;
      p | m_rangeFieldFloor;
      p | m_rangeHistFloor;
      p | m_rangeIntegFloor;
      p | m_dt;
      p | m_dtn;
      p | m_nvol;
      p | m_boxnodes;
      p | m_disc;
      p | m_transporter;
      p | m_meshwriter;
      p | m_refiner;
      p | m_el;
      if (p.isUnpacking()) {
        m_inpoel = std::get< 0 >( m_el );
        m_gid = std::get< 1 >( m_el );
        m_lid = std::get< 2 >( m_el );
      }
      p | m_coord;
      p | m_nodeCommMap;
      p | m_meshvol;
      p | m_v;
      p | m_vol;
      p | m_volc;
      p | m_cvolc;
      p | m_boxvol;
      p | m_timer;
      p | m_refined;
      p( reinterpret_cast<char*>(&m_prevstatus), sizeof(Clock::time_point) );
      p | m_nrestart;
      p | m_histdata;
      p | m_res;
      p | m_res0;
      p | m_res1;
      p | m_dea;
      p | m_deastarted;
      p | m_pit;
      p | m_mit;
    }
    //! \brief Pack/Unpack serialize operator|
    //! \param[in,out] p Charm++'s PUP::er serializer object reference
    //! \param[in,out] i Discretization object reference
    friend void operator|( PUP::er& p, Discretization& i ) { i.pup(p); }
    //@}

  private:
    // Shorthand for clock, setting an internal clock type
    using Clock = std::chrono::high_resolution_clock;

    //! Total number of mesh points (across all meshes)
    std::size_t m_npoin;
    //! Mesh ID
    std::size_t m_meshid;
    //! Function to call after mesh-to-mesh solution transfer is complete
    CkCallback m_transfer_complete;
    //! Pointer to solution during mesh-to-mesh solution transfer
    tk::Fields* m_transfer_sol;
    //! Transfer flags at nodes (if coupled)
    std::vector< double > m_transfer_flag;
    //! Transfer flags if true
    bool m_trflag;
    //! Holes tessellation
    //! \details: Key: hole id, value: list of triangle node coordinates
    std::unordered_map< std::size_t, std::vector< tk::real > > m_transfer_hole;
    //! Callback to continue with after hole communication
    CkCallback m_holcont;
    //! Total number of Discretization chares
    int m_nchare;
    //! Iteration count
    uint64_t m_it;
    //! Iteration count with mesh refinement
    //! \details Used as the restart sequence number {RS} in saving output in
    //!    an ExodusII sequence
    //! \see https://www.paraview.org/Wiki/Restarted_Simulation_Readers
    uint64_t m_itr;
    //! Field output iteration count without mesh refinement
    //! \details Counts the number of field outputs to file during two
    //!   time steps with mesh efinement
    uint64_t m_itf;
    //! Flag that is nonzero during setup and zero during time stepping
    std::size_t m_initial;
    //! Physical time
    tk::real m_t;
    //! Physics time at last field output
    tk::real m_lastDumpTime;
    //! Recent floor of physics time divided by field output interval time
    tk::real m_physFieldFloor;
    //! Recent floor of physics time divided by history output interval time
    tk::real m_physHistFloor;
    //! Recent floor of physics time divided by integral output interval time
    tk::real m_physIntegFloor;
    //! Recent floors of physics time divided by field output time for ranges
    std::vector< tk::real > m_rangeFieldFloor;
    //! Recent floors of physics time divided by history output time for ranges
    std::vector< tk::real > m_rangeHistFloor;
    //! Recent floors of physics time divided by integral output time for ranges
    std::vector< tk::real > m_rangeIntegFloor;
    //! Physical time step size
    tk::real m_dt;
    //! Physical time step size at the previous time step
    tk::real m_dtn;
    //! \brief Number of chares from which we received nodal volume
    //!   contributions on chare boundaries
    std::size_t m_nvol;
    //! List of nodes at which box user ICs are set for each IC box
    std::vector< std::unordered_set< std::size_t > > m_boxnodes;
    //! Discretization proxy for all meshes
    std::vector< CProxy_Discretization > m_disc;
    //! Transporter proxy
    CProxy_Transporter m_transporter;
    //! Mesh writer proxy
    tk::CProxy_MeshWriter m_meshwriter;
    //! Mesh refiner proxy
    CProxy_Refiner m_refiner;
    //! \brief Elements of the mesh chunk we operate on
    //! \details Initialized by the constructor. The first vector is the element
    //!   connectivity (local IDs), the second vector is the global node IDs of
    //!   owned elements, while the third one is a map of global->local node
    //!   IDs.
    tk::UnsMesh::Chunk m_el;
    //! Alias to element connectivity
    std::vector< std::size_t >& m_inpoel = std::get<0>( m_el );
    //! Alias to global node IDs of owned elements
    std::vector< std::size_t >& m_gid = std::get<1>( m_el );
    //! \brief Alias to local node ids associated to the global ones of owned
    //!    elements
    std::unordered_map< std::size_t, std::size_t >& m_lid = std::get<2>( m_el );
    //! Mesh point coordinates
    tk::UnsMesh::Coords m_coord;
    //! \brief Global mesh node IDs bordering the mesh chunk held by fellow
    //!   Discretization chares associated to their chare IDs
    std::unordered_map< int, std::unordered_set< std::size_t > > m_nodeCommMap;
    //! Total mesh volume
    tk::real m_meshvol;
    //! Nodal mesh volumes
    //! \details This is the volume of the mesh associated to nodes of owned
    //!   elements (sum of surrounding cell volumes / 4) without contributions
    //!   from other chares on chare-boundaries
    std::vector< tk::real > m_v;
    //! Volume of nodes
    //! \details This is the volume of the mesh associated to nodes of owned
    //!   elements (sum of surrounding cell volumes / 4) with contributions from
    //!   other chares on chare-boundaries
    std::vector< tk::real > m_vol;
    //! Receive buffer for volume of nodes (with global node id as key)
    //! \details This is a communication buffer used to compute the volume of
    //!   the mesh associated to nodes of owned elements (sum of surrounding
    //!   cell volumes / 4) with contributions from other chares on
    //!   chare-boundaries.
    std::unordered_map< std::size_t, tk::real > m_volc;
    //! Chare-volume contributions per chare per chare-boundary node
    //! \details Outer key: chare id, value: local node id and node volume
    std::unordered_map< int, std::unordered_map< std::size_t, tk::real > >
      m_cvolc;
    //! Volume of user-defined box IC
    tk::real m_boxvol;
    //! Timer measuring a time step
    tk::Timer m_timer;
    //! 1 if mesh was refined in a time step, 0 if it was not
    int m_refined;
    //! Time point storing clock state at status()
    Clock::time_point m_prevstatus;
    //! Number of times restarted
    int m_nrestart;
    //! Data at history point locations
    std::vector< HistData > m_histdata;
    //! Current residual (during convergence to steady state)
    tk::real m_res;
    //! Residual at previous ETA calculation (during convergence to steady state)
    tk::real m_res0;
    //! Residual at next ETA calculation (during convergence to steady state)
    tk::real m_res1;
    //! Numberf of deactivated chares
    int m_dea;
    //! Flag: 1 if deactivation procedure has started
    int m_deastarted;
    //! Numberf of recent pressure solve iterations taken
    std::size_t m_pit;
    //! Numberf of recent momentum/transport solve iterations taken
    std::size_t m_mit;

    //! Set mesh coordinates based on coordinates map
    tk::UnsMesh::Coords setCoord( const tk::UnsMesh::CoordMap& coordmap );
};

} // inciter::