Xyst test code coverage report
Current view: top level - Inciter - Transporter.cpp (source / functions) Hit Total Coverage
Commit: d790e211db0f2cf155e72db869cf1a5a372c10f5 Lines: 703 770 91.3 %
Date: 2025-02-17 15:57:47 Functions: 41 44 93.2 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 824 1784 46.2 %

           Branch data     Line data    Source code
       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/Inciter/Transporter.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                 :            :   \brief     Transporter drives the time integration of transport equations
      10                 :            :   \details   Transporter drives the time integration of transport equations.
      11                 :            :     The implementation uses the Charm++ runtime system and is fully asynchronous,
      12                 :            :     overlapping computation, communication as well as I/O. The algorithm
      13                 :            :     utilizes the structured dagger (SDAG) Charm++ functionality. The high-level
      14                 :            :     overview of the algorithm structure and how it interfaces with Charm++ is
      15                 :            :     discussed in the Charm++ interface file src/Inciter/transporter.ci.
      16                 :            : */
      17                 :            : // *****************************************************************************
      18                 :            : 
      19                 :            : #include <string>
      20                 :            : #include <iomanip>
      21                 :            : #include <cstddef>
      22                 :            : #include <unordered_set>
      23                 :            : 
      24                 :            : #include "Transporter.hpp"
      25                 :            : #include "Fields.hpp"
      26                 :            : #include "UniPDF.hpp"
      27                 :            : #include "PDFWriter.hpp"
      28                 :            : #include "ContainerUtil.hpp"
      29                 :            : #include "LoadDistributor.hpp"
      30                 :            : #include "ExodusIIMeshReader.hpp"
      31                 :            : #include "InciterConfig.hpp"
      32                 :            : #include "DiagWriter.hpp"
      33                 :            : #include "Diagnostics.hpp"
      34                 :            : #include "Integrals.hpp"
      35                 :            : #include "Callback.hpp"
      36                 :            : #include "Problems.hpp"
      37                 :            : 
      38                 :            : #include "NoWarning/inciter.decl.h"
      39                 :            : #include "NoWarning/partitioner.decl.h"
      40                 :            : 
      41                 :            : extern CProxy_Main mainProxy;
      42                 :            : 
      43                 :            : namespace inciter {
      44                 :            : 
      45                 :            : extern ctr::Config g_cfg;
      46                 :            : extern int g_nrestart;
      47                 :            : 
      48                 :            : }
      49                 :            : 
      50                 :            : using inciter::Transporter;
      51                 :            : 
      52                 :        254 : Transporter::Transporter() :
      53 [ +  - ][ +  + ]:        508 :   m_input{ g_cfg.get< tag::input >() },
                 [ -  - ]
      54         [ +  - ]:        254 :   m_nchare( m_input.size() ),
      55         [ +  - ]:        254 :   m_ncit( m_nchare.size(), 0 ),
      56                 :        254 :   m_nload( 0 ),
      57                 :        254 :   m_npart( 0 ),
      58                 :        254 :   m_nstat( 0 ),
      59                 :        254 :   m_ndisc( 0 ),
      60                 :        254 :   m_nchk( 0 ),
      61                 :        254 :   m_ncom( 0 ),
      62         [ +  - ]:        254 :   m_nt0refit( m_nchare.size(), 0 ),
      63         [ +  - ]:        254 :   m_ndtrefit( m_nchare.size(), 0 ),
      64         [ +  - ]:        254 :   m_noutrefit( m_nchare.size(), 0 ),
      65         [ +  - ]:        254 :   m_noutderefit( m_nchare.size(), 0 ),
      66         [ +  - ]:        254 :   m_nelem( m_nchare.size() ),
      67         [ +  - ]:        254 :   m_finished( m_nchare.size(), 0 ),
      68         [ +  - ]:        254 :   m_meshvol( m_nchare.size() ),
      69         [ +  - ]:        254 :   m_minstat( m_nchare.size() ),
      70         [ +  - ]:        254 :   m_maxstat( m_nchare.size() ),
      71         [ +  - ]:        254 :   m_avgstat( m_nchare.size() ),
      72         [ +  - ]:        254 :   m_progMesh( g_cfg.get< tag::feedback >(), ProgMeshPrefix, ProgMeshLegend ),
      73         [ +  - ]:       1524 :   m_progWork( g_cfg.get< tag::feedback >(), ProgWorkPrefix, ProgWorkLegend )
      74                 :            : // *****************************************************************************
      75                 :            : //  Constructor
      76                 :            : // *****************************************************************************
      77                 :            : {
      78                 :        254 :   const auto nstep = g_cfg.get< tag::nstep >();
      79                 :        254 :   const auto t0 = g_cfg.get< tag::t0 >();
      80                 :        254 :   const auto term = g_cfg.get< tag::term >();
      81                 :        254 :   const auto constdt = g_cfg.get< tag::dt >();
      82                 :            : 
      83                 :            :   // If the desired max number of time steps is larger than zero, and the
      84                 :            :   // termination time is larger than the initial time, and the constant time
      85                 :            :   // step size (if that is used) is smaller than the duration of the time to be
      86                 :            :   // simulated, we have work to do, otherwise, finish right away. If a constant
      87                 :            :   // dt is not used, that part of the logic is always true as the default
      88                 :            :   // constdt is zero.
      89 [ +  - ][ +  - ]:        254 :   if ( nstep != 0 && term > t0 && constdt < term-t0 ) {
                 [ +  - ]
      90                 :            : 
      91                 :            :     // Enable SDAG waits for collecting mesh statistics
      92         [ +  - ]:        254 :     thisProxy.wait4stat();
      93                 :            : 
      94                 :            :     // Configure and write diagnostics file header
      95         [ +  - ]:        254 :     diagHeader();
      96                 :            : 
      97                 :            :     // Configure and write integrals file header
      98         [ +  - ]:        254 :     integralsHeader();
      99                 :            : 
     100                 :            :     // Create mesh partitioner AND boundary condition object group
     101         [ +  - ]:        254 :     createPartitioner();
     102                 :            : 
     103         [ -  - ]:          0 :   } else finish();      // stop if no time stepping requested
     104 [ +  - ][ -  - ]:        507 : }
                 [ -  - ]
     105                 :            : 
     106                 :         12 : Transporter::Transporter( CkMigrateMessage* m ) :
     107                 :            :   CBase_Transporter( m ),
     108         [ +  - ]:         12 :   m_progMesh( g_cfg.get< tag::feedback >(), ProgMeshPrefix, ProgMeshLegend ),
     109         [ +  - ]:         24 :   m_progWork( g_cfg.get< tag::feedback >(), ProgWorkPrefix, ProgWorkLegend )
     110                 :            : // *****************************************************************************
     111                 :            : //  Migrate constructor: returning from a checkpoint
     112                 :            : //! \param[in] m Charm++ migrate message
     113                 :            : // *****************************************************************************
     114                 :            : {
     115                 :         12 :   auto print = tk::Print();
     116         [ +  - ]:         12 :   print << "\nXyst> Restarted from checkpoint\n";
     117         [ +  - ]:         12 :   inthead( print );
     118                 :         12 : }
     119                 :            : 
     120                 :            : bool
     121                 :        253 : Transporter::matchBCs( std::map< int, std::vector< std::size_t > >& bnd )
     122                 :            : // *****************************************************************************
     123                 :            : // Verify that side sets specified in the control file exist in mesh file
     124                 :            : //! \details This function does two things: (1) it verifies that the side
     125                 :            : //!   sets used in the input file (either to which boundary conditions (BC)
     126                 :            : //!   are assigned or listed as field output by the user in the
     127                 :            : //!   input file) all exist among the side sets read from the input mesh
     128                 :            : //!   file and errors out if at least one does not, and (2) it matches the
     129                 :            : //!   side set ids at which the user has configured BCs (or listed as an output
     130                 :            : //!   surface) to side set ids read from the mesh file and removes those face
     131                 :            : //!   and node lists associated to side sets that the user did not set BCs or
     132                 :            : //!   listed as field output on (as they will not need processing further since
     133                 :            : //!   they will not be used).
     134                 :            : //! \param[in,out] bnd Node or face lists mapped to side set ids
     135                 :            : //! \return True if sidesets have been used and found in mesh
     136                 :            : // *****************************************************************************
     137                 :            : {
     138                 :        253 :   std::unordered_set< int > usersets;
     139                 :            : 
     140                 :            :   // Collect side sets at which BCs are set
     141                 :            : 
     142         [ +  + ]:        669 :   for (const auto& s : g_cfg.get< tag::bc_dir >()) {
     143 [ +  - ][ +  - ]:        416 :     if (!s.empty()) usersets.insert( s[0] );
     144                 :            :   }
     145                 :            : 
     146 [ +  - ][ +  + ]:        506 :   for (auto s : g_cfg.get< tag::bc_sym >()) usersets.insert( s );
     147                 :            : 
     148 [ +  - ][ +  + ]:        263 :   for (auto s : g_cfg.get< tag::bc_far >()) usersets.insert( s );
     149                 :            :  
     150         [ +  + ]:        265 :   for (const auto& s : g_cfg.get< tag::bc_pre >()) {
     151 [ +  - ][ +  - ]:         12 :     if (!s.empty()) usersets.insert( s[0] );
     152                 :            :   }
     153                 :            :  
     154         [ +  + ]:        325 :   for (const auto& s : g_cfg.get< tag::pre_bc_dir >()) {
     155 [ +  - ][ +  - ]:         72 :     if (!s.empty()) usersets.insert( s[0] );
     156                 :            :   }
     157                 :            : 
     158 [ +  - ][ +  + ]:        256 :   for (auto s : g_cfg.get< tag::pre_bc_sym >()) usersets.insert( s );
     159                 :            : 
     160                 :            :   // Add sidesets requested for field output
     161 [ +  - ][ +  + ]:        286 :   for (auto s : g_cfg.get< tag::fieldout >()) usersets.insert( s );
     162                 :            :   // Add sidesets requested for integral output
     163 [ +  - ][ +  + ]:        266 :   for (auto s : g_cfg.get< tag::integout >()) usersets.insert( s );
     164                 :            : 
     165                 :            :   // Find user-configured side set ids among side sets read from mesh file
     166                 :        253 :   std::unordered_set< int > sidesets_used;
     167         [ +  + ]:       1012 :   for (auto i : usersets) {       // for all side sets used in control file
     168 [ +  - ][ +  - ]:        759 :     if (bnd.find(i) != end(bnd))  // used set found among side sets in file
     169         [ +  - ]:        759 :       sidesets_used.insert( i );  // store side set id configured as BC
     170                 :            :     else {
     171 [ -  - ][ -  - ]:          0 :       Throw( "Boundary conditions specified on side set " + std::to_string(i) +
         [ -  - ][ -  - ]
     172                 :            :              " which does not exist in mesh file" );
     173                 :            :     }
     174                 :            :   }
     175                 :            :  
     176                 :            :   // Remove sidesets not used (will not process those further)
     177         [ +  - ]:        253 :   tk::erase_if( bnd, [&]( auto& item ) {
     178         [ +  - ]:       1382 :     return sidesets_used.find( item.first ) == end(sidesets_used);
     179                 :            :   });
     180                 :            : 
     181                 :        506 :   return not bnd.empty();
     182                 :        253 : }
     183                 :            : 
     184                 :            : void
     185                 :        254 : Transporter::createPartitioner()
     186                 :            : // *****************************************************************************
     187                 :            : // Create mesh partitioner AND boundary conditions group
     188                 :            : // *****************************************************************************
     189                 :            : {
     190                 :            :   // cppcheck-suppress unreadVariable
     191                 :        254 :   auto print = tk::Print();
     192                 :            : 
     193                 :            :   // Create partitioner callbacks (order important)
     194                 :            :   tk::PartitionerCallback cbp {{
     195 [ +  - ][ +  - ]:        254 :       CkCallback( CkReductionTarget(Transporter,queriedPart), thisProxy )
     196 [ +  - ][ +  - ]:        508 :     , CkCallback( CkReductionTarget(Transporter,respondedPart), thisProxy )
     197 [ +  - ][ +  - ]:        508 :     , CkCallback( CkReductionTarget(Transporter,load), thisProxy )
     198 [ +  - ][ +  - ]:        508 :     , CkCallback( CkReductionTarget(Transporter,partitioned), thisProxy )
     199 [ +  - ][ +  - ]:        508 :     , CkCallback( CkReductionTarget(Transporter,distributed), thisProxy )
     200 [ +  - ][ +  - ]:        508 :     , CkCallback( CkReductionTarget(Transporter,refinserted), thisProxy )
     201                 :        254 :   }};
     202                 :            : 
     203                 :            :   // Create refiner callbacks (order important)
     204                 :            :   tk::RefinerCallback cbr {{
     205 [ +  - ][ +  - ]:        254 :       CkCallback( CkReductionTarget(Transporter,queriedRef), thisProxy )
     206 [ +  - ][ +  - ]:        508 :     , CkCallback( CkReductionTarget(Transporter,respondedRef), thisProxy )
     207 [ +  - ][ +  - ]:        508 :     , CkCallback( CkReductionTarget(Transporter,compatibility), thisProxy )
     208 [ +  - ][ +  - ]:        508 :     , CkCallback( CkReductionTarget(Transporter,bndint), thisProxy )
     209 [ +  - ][ +  - ]:        508 :     , CkCallback( CkReductionTarget(Transporter,matched), thisProxy )
     210 [ +  - ][ +  - ]:        508 :     , CkCallback( CkReductionTarget(Transporter,refined), thisProxy )
     211                 :        254 :   }};
     212                 :            : 
     213                 :            :   // Create sorter callbacks (order important)
     214                 :            :   tk::SorterCallback cbs {{
     215 [ +  - ][ +  - ]:        254 :       CkCallback( CkReductionTarget(Transporter,queried), thisProxy )
     216 [ +  - ][ +  - ]:        508 :     , CkCallback( CkReductionTarget(Transporter,responded), thisProxy )
     217 [ +  - ][ +  - ]:        508 :     , CkCallback( CkReductionTarget(Transporter,discinserted), thisProxy )
     218 [ +  - ][ +  - ]:        508 :     , CkCallback( CkReductionTarget(Transporter,workinserted), thisProxy )
     219                 :        254 :   }};
     220                 :            : 
     221                 :            :   // Start timer measuring preparation of mesh(es) for partitioning
     222         [ +  - ]:        254 :   m_timer[ TimerTag::MESH_READ ];
     223                 :            : 
     224 [ -  + ][ -  - ]:        254 :   ErrChk( !m_input.empty(), "No input mesh" );
         [ -  - ][ -  - ]
     225                 :            : 
     226                 :            :   // Start preparing mesh(es)
     227 [ -  + ][ +  - ]:        254 :   print.section( "Reading mesh" + std::string(m_input.size()>1?"es":"") );
         [ +  - ][ +  - ]
     228                 :            : 
     229                 :            :   // Read boundary (side set) data from a list of input mesh files
     230                 :        254 :   std::size_t meshid = 0;
     231         [ +  + ]:        507 :   for (const auto& filename : m_input) {
     232                 :            :     // Create mesh reader for reading side sets from file
     233         [ +  - ]:        254 :     tk::ExodusIIMeshReader mr( filename );
     234                 :            : 
     235                 :            :     // Read out total number of mesh points from mesh file
     236 [ +  - ][ +  - ]:        254 :     m_npoin.push_back( mr.npoin() );
     237                 :            : 
     238                 :        254 :     std::map< int, std::vector< std::size_t > > bface;
     239                 :        254 :     std::map< int, std::vector< std::size_t > > faces;
     240                 :        254 :     std::map< int, std::vector< std::size_t > > bnode;
     241                 :            : 
     242                 :            :     // Read boundary-face connectivity on side sets
     243         [ +  - ]:        254 :     mr.readSidesetFaces( bface, faces );
     244                 :            : 
     245                 :        253 :     bool bcs_set = false;
     246                 :            :     // Read node lists on side sets
     247         [ +  - ]:        253 :     bnode = mr.readSidesetNodes();
     248                 :            :     // Verify boundarty condition (BC) side sets used exist in mesh file
     249         [ +  - ]:        253 :     bcs_set = matchBCs( bnode );
     250 [ -  + ][ -  - ]:        253 :     bcs_set = bcs_set || matchBCs( bface );
                 [ -  - ]
     251                 :            : 
     252                 :            :     // Warn on no BCs
     253 [ -  + ][ -  - ]:        253 :     if (!bcs_set) print << "\n>>> WARNING: No boundary conditions set\n\n";
     254                 :            : 
     255                 :            :     // Create empty discretization chare array
     256 [ +  - ][ +  - ]:        253 :     m_discretization.push_back( CProxy_Discretization::ckNew() );
                 [ +  - ]
     257         [ +  - ]:        253 :     CkArrayOptions opt;
     258 [ +  - ][ +  - ]:        253 :     opt.bindTo( m_discretization.back() );
     259                 :            : 
     260                 :            :     // Create empty discretization scheme chare array (bound to discretization)
     261                 :        253 :     CProxy_RieCG riecg;
     262                 :        253 :     CProxy_LaxCG laxcg;
     263                 :        253 :     CProxy_ZalCG zalcg;
     264                 :        253 :     CProxy_KozCG kozcg;
     265                 :        253 :     CProxy_ChoCG chocg;
     266                 :        253 :     CProxy_LohCG lohcg;
     267                 :        253 :     tk::CProxy_ConjugateGradients cgpre, cgmom;
     268                 :        253 :     const auto& solver = g_cfg.get< tag::solver >();
     269         [ +  + ]:        253 :     if (solver == "riecg") {
     270 [ +  - ][ +  - ]:         73 :       m_riecg.push_back( CProxy_RieCG::ckNew(opt) );
                 [ +  - ]
     271         [ +  - ]:         73 :       riecg = m_riecg.back();
     272                 :            :     }
     273         [ +  + ]:        180 :     else if (solver == "laxcg") {
     274 [ +  - ][ +  - ]:         23 :       m_laxcg.push_back( CProxy_LaxCG::ckNew(opt) );
                 [ +  - ]
     275         [ +  - ]:         23 :       laxcg = m_laxcg.back();
     276                 :            :     }
     277         [ +  + ]:        157 :     else if (solver == "zalcg") {
     278 [ +  - ][ +  - ]:         36 :       m_zalcg.push_back( CProxy_ZalCG::ckNew(opt) );
                 [ +  - ]
     279         [ +  - ]:         36 :       zalcg = m_zalcg.back();
     280                 :            :     }
     281         [ +  + ]:        121 :     else if (solver == "kozcg") {
     282 [ +  - ][ +  - ]:         44 :       m_kozcg.push_back( CProxy_KozCG::ckNew(opt) );
                 [ +  - ]
     283         [ +  - ]:         44 :       kozcg = m_kozcg.back();
     284                 :            :     }
     285         [ +  + ]:         77 :     else if (solver == "chocg") {
     286 [ +  - ][ +  - ]:         44 :       m_chocg.push_back( CProxy_ChoCG::ckNew(opt) );
                 [ +  - ]
     287         [ +  - ]:         44 :       chocg = m_chocg.back();
     288 [ +  - ][ +  - ]:         44 :       m_cgpre.push_back( tk::CProxy_ConjugateGradients::ckNew(opt) );
                 [ +  - ]
     289         [ +  - ]:         44 :       cgpre = m_cgpre.back();
     290 [ +  - ][ +  - ]:         44 :       m_cgmom.push_back( tk::CProxy_ConjugateGradients::ckNew(opt) );
                 [ +  - ]
     291         [ +  - ]:         44 :       cgmom = m_cgmom.back();
     292                 :            :     }
     293         [ +  - ]:         33 :     else if (solver == "lohcg") {
     294 [ +  - ][ +  - ]:         33 :       m_lohcg.push_back( CProxy_LohCG::ckNew(opt) );
                 [ +  - ]
     295         [ +  - ]:         33 :       lohcg = m_lohcg.back();
     296 [ +  - ][ +  - ]:         33 :       m_cgpre.push_back( tk::CProxy_ConjugateGradients::ckNew(opt) );
                 [ +  - ]
     297         [ +  - ]:         33 :       cgpre = m_cgpre.back();
     298                 :            :     }
     299                 :            :     else {
     300 [ -  - ][ -  - ]:          0 :       Throw( "Unknown solver: " + solver );
                 [ -  - ]
     301                 :            :     }
     302                 :            : 
     303                 :            :     // Create empty mesh refiner chare array (bound to discretization)
     304 [ +  - ][ +  - ]:        253 :     m_refiner.push_back( CProxy_Refiner::ckNew(opt) );
                 [ +  - ]
     305                 :            :     // Create empty mesh sorter Charm++ chare array (bound to discretization)
     306 [ +  - ][ +  - ]:        253 :     m_sorter.push_back( CProxy_Sorter::ckNew(opt) );
                 [ +  - ]
     307                 :            : 
     308                 :            :     // Create MeshWriter chare group for mesh
     309 [ +  - ][ +  - ]:        253 :     m_meshwriter.push_back(
     310                 :            :       tk::CProxy_MeshWriter::ckNew(
     311                 :        253 :         g_cfg.get< tag::benchmark >(), m_input.size() ) );
     312                 :            : 
     313                 :            :     // Create mesh partitioner Charm++ chare nodegroup for all meshes
     314         [ +  - ]:        253 :     m_partitioner.push_back(
     315                 :            :       CProxy_Partitioner::ckNew( meshid, filename, cbp, cbr, cbs,
     316         [ +  - ]:        253 :         thisProxy, m_refiner.back(), m_sorter.back(), m_meshwriter.back(),
     317                 :        253 :         m_discretization.back(), riecg, laxcg, zalcg, kozcg, chocg, lohcg,
     318                 :            :         cgpre, cgmom, bface, faces, bnode ) );
     319                 :            : 
     320                 :        253 :     ++meshid;
     321                 :        253 :   }
     322                 :        253 : }
     323                 :            : 
     324                 :            : void
     325                 :        253 : Transporter::load( std::size_t meshid, std::size_t nelem )
     326                 :            : // *****************************************************************************
     327                 :            : // Reduction target: the mesh has been read from file on all PEs
     328                 :            : //! \param[in] meshid Mesh id (summed accross all compute nodes)
     329                 :            : //! \param[in] nelem Number of mesh elements per mesh (summed across all
     330                 :            : //!    compute nodes)
     331                 :            : // *****************************************************************************
     332                 :            : {
     333                 :        253 :   meshid /= static_cast< std::size_t >( CkNumNodes() );
     334 [ -  + ][ -  - ]:        253 :   Assert( meshid < m_nelem.size(), "MeshId indexing out" );
         [ -  - ][ -  - ]
     335                 :        253 :   m_nelem[meshid] = nelem;
     336                 :            : 
     337                 :            :   // Compute load distribution given total work (nelem) and user-specified
     338                 :            :   // virtualization
     339                 :            :   uint64_t chunksize, remainder;
     340                 :        506 :   m_nchare[meshid] = static_cast<int>(
     341         [ +  - ]:        506 :     tk::linearLoadDistributor(
     342                 :        253 :        g_cfg.get< tag::virt >(),
     343                 :        253 :        m_nelem[meshid], CkNumPes(), chunksize, remainder ) );
     344                 :            : 
     345                 :            :   // Store sum of meshids (across all chares, key) for each meshid (value).
     346                 :            :   // This is used to look up the mesh id after collectives that sum their data.
     347         [ +  - ]:        253 :   m_meshid[ static_cast<std::size_t>(m_nchare[meshid])*meshid ] = meshid;
     348 [ -  + ][ -  - ]:        253 :   Assert( meshid < m_nelem.size(), "MeshId indexing out" );
         [ -  - ][ -  - ]
     349                 :            : 
     350                 :            :   // Partition first mesh
     351         [ +  - ]:        253 :   if (meshid == 0) {
     352         [ +  - ]:        253 :     m_timer[ TimerTag::MESH_PART ];  // start timer measuring mesh partitioning
     353         [ +  - ]:        253 :     m_partitioner[0].partition( m_nchare[0] );
     354                 :            :   }
     355                 :            : 
     356         [ +  - ]:        253 :   if (++m_nload == m_nelem.size()) {     // all meshes have been loaded
     357                 :        253 :     m_nload = 0;
     358                 :        253 :     auto print = tk::Print();
     359                 :            : 
     360         [ +  - ]:        253 :     auto& timer = tk::ref_find( m_timer, TimerTag::MESH_READ );
     361         [ +  - ]:        253 :     timer.second = timer.first.dsec();
     362 [ +  - ][ +  - ]:        253 :     print << "Mesh read time: " + std::to_string( timer.second ) + " sec\n";
         [ +  - ][ +  - ]
     363                 :            : 
     364                 :            :     // Print out mesh partitioning configuration
     365 [ +  - ][ +  - ]:        253 :     print.section( "Partitioning mesh" );
     366 [ +  - ][ +  - ]:        253 :     print.item( "Partitioner", g_cfg.get< tag::part >() );
     367 [ +  - ][ +  - ]:        253 :     print.item( "Virtualization", g_cfg.get< tag::virt >() );
     368                 :            :     // Print out initial mesh statistics
     369 [ +  - ][ +  - ]:        253 :     meshstat( "Mesh read from file" );
     370                 :            : 
     371                 :            :     // Tell meshwriter the total number of chares
     372         [ +  - ]:        253 :     m_meshwriter[meshid].nchare( m_nchare[meshid] );
     373                 :            : 
     374                 :            :     // Query number of initial mesh refinement steps
     375                 :        253 :     int nref = 0;
     376         [ +  + ]:        253 :     if (g_cfg.get< tag::href_t0 >()) {
     377                 :         12 :       nref = static_cast<int>(g_cfg.get< tag::href_init >().size());
     378                 :            :     }
     379                 :            : 
     380                 :            :     // Query if PE-local reorder is configured
     381                 :        253 :     int nreord = 0;
     382         [ +  + ]:        253 :     if (g_cfg.get< tag::reorder >()) nreord = m_nchare[0];
     383                 :            : 
     384         [ +  - ]:        253 :     print << '\n';
     385 [ +  - ][ +  - ]:        506 :     m_progMesh.start( print, "Preparing mesh", {{ CkNumPes(), CkNumPes(), nref,
     386                 :        253 :       m_nchare[0], m_nchare[0], nreord, nreord }} );
     387                 :            :   }
     388                 :        253 : }
     389                 :            : 
     390                 :            : void
     391                 :        253 : Transporter::partitioned( std::size_t meshid )
     392                 :            : // *****************************************************************************
     393                 :            : // Reduction target: a mesh has been partitioned
     394                 :            : //! \param[in] meshid Mesh id
     395                 :            : // *****************************************************************************
     396                 :            : {
     397         [ +  - ]:        253 :   if (++m_npart == m_nelem.size()) {     // all meshes have been partitioned
     398                 :        253 :     m_npart = 0;
     399         [ +  - ]:        253 :     auto& timer = tk::ref_find( m_timer, TimerTag::MESH_PART );
     400                 :        253 :     timer.second = timer.first.dsec();
     401         [ +  - ]:        253 :     m_timer[ TimerTag::MESH_DIST ];  // start timer measuring mesh distribution
     402                 :            :   } else { // partition next mesh
     403                 :          0 :     m_partitioner[meshid+1].partition( m_nchare[meshid+1] );
     404                 :            :   }
     405                 :        253 : }
     406                 :            : 
     407                 :            : void
     408                 :        253 : Transporter::distributed( std::size_t meshid )
     409                 :            : // *****************************************************************************
     410                 :            : // Reduction target: all compute nodes have distributed their mesh after
     411                 :            : // partitioning
     412                 :            : //! \param[in] meshid Mesh id
     413                 :            : // *****************************************************************************
     414                 :            : {
     415                 :        253 :   m_partitioner[meshid].refine();
     416         [ +  - ]:        253 :   auto& timer = tk::ref_find( m_timer, TimerTag::MESH_DIST );
     417                 :        253 :   timer.second = timer.first.dsec();
     418                 :        253 : }
     419                 :            : 
     420                 :            : void
     421                 :        253 : Transporter::refinserted( std::size_t meshid, std::size_t error )
     422                 :            : // *****************************************************************************
     423                 :            : // Reduction target: all compute nodes have created the mesh refiners
     424                 :            : //! \param[in] meshid Mesh id (aggregated across all compute nodes with operator
     425                 :            : //!   max)
     426                 :            : //! \param[in] error Error code (aggregated across all compute nodes with
     427                 :            : //!   operator max)
     428                 :            : // *****************************************************************************
     429                 :            : {
     430         [ -  + ]:        253 :   if (error) {
     431                 :            : 
     432                 :          0 :     tk::Print() <<
     433                 :            :         "\n>>> ERROR: A worker chare was not assigned any mesh "
     434 [ -  - ][ -  - ]:          0 :         "elements after distributing mesh " + std::to_string(meshid) +
                 [ -  - ]
     435                 :            :         ". This can happen in SMP-mode with a large +ppn parameter (number "
     436                 :            :         "of worker threads per logical node) or using zoltan's hypergraph "
     437                 :            :         "partitioning (phg), which is non-determinsitic. Solution 1: In SMP "
     438                 :            :         "mode decrease +ppn. Solution 2: Try a different partitioning "
     439                 :            :         "algorithm, e.g., rcb, rib, or hsfc, or configure phg differently by "
     440                 :            :         "passing extra zoltan parameters in the control file. To learn how, "
     441         [ -  - ]:          0 :         "grep tests/ for 'zoltan_params'.";
     442                 :          0 :     finish( meshid );
     443                 :            : 
     444                 :            :   } else {
     445                 :            : 
     446                 :        253 :     m_refiner[meshid].doneInserting();
     447                 :            : 
     448                 :            :   }
     449                 :        253 : }
     450                 :            : 
     451                 :            : void
     452                 :         38 : Transporter::queriedRef( std::size_t meshid )
     453                 :            : // *****************************************************************************
     454                 :            : // Reduction target: all Refiner chares have queried their boundary edges
     455                 :            : //! \param[in] meshid Mesh id
     456                 :            : // *****************************************************************************
     457                 :            : {
     458                 :         38 :   m_refiner[meshid].response();
     459                 :         38 : }
     460                 :            : 
     461                 :            : void
     462                 :         38 : Transporter::respondedRef( std::size_t meshid )
     463                 :            : // *****************************************************************************
     464                 :            : // Reduction target: all Refiner chares have setup their boundary edges
     465                 :            : //! \param[in] meshid Mesh id
     466                 :            : // *****************************************************************************
     467                 :            : {
     468                 :         38 :   m_refiner[meshid].refine();
     469                 :         38 : }
     470                 :            : 
     471                 :            : void
     472                 :         25 : Transporter::compatibility( std::size_t meshid )
     473                 :            : // *****************************************************************************
     474                 :            : // Reduction target: all Refiner chares have received a round of edges,
     475                 :            : // and have run their compatibility algorithm
     476                 :            : //! \param[in] meshid Mesh id (aggregated across all chares using operator max)
     477                 :            : //! \details This is called iteratively, until convergence by Refiner. At this
     478                 :            : //!   point all Refiner chares have received a round of edge data (tags whether
     479                 :            : //!   an edge needs to be refined, etc.), and applied the compatibility
     480                 :            : //!   algorithm independent of other Refiner chares. We keep going until the
     481                 :            : //!   mesh is no longer modified by the compatibility algorithm, based on a new
     482                 :            : //!   round of edge data communication started in Refiner::comExtra().
     483                 :            : // *****************************************************************************
     484                 :            : {
     485                 :         25 :   m_refiner[meshid].correctref();
     486                 :         25 : }
     487                 :            : 
     488                 :            : void
     489                 :         38 : Transporter::matched( std::size_t summeshid,
     490                 :            :                       std::size_t nextra,
     491                 :            :                       std::size_t nref,
     492                 :            :                       std::size_t nderef,
     493                 :            :                       std::size_t sumrefmode )
     494                 :            : // *****************************************************************************
     495                 :            : // Reduction target: all Refiner chares have matched/corrected the tagging
     496                 :            : // of chare-boundary edges, all chares are ready to perform refinement.
     497                 :            : //! \param[in] summeshid Mesh id (summed across all chares)
     498                 :            : //! \param[in] nextra Sum (across all chares) of the number of edges on each
     499                 :            : //!   chare that need correction along chare boundaries
     500                 :            : //! \param[in] nref Sum of number of refined tetrahedra across all chares.
     501                 :            : //! \param[in] nderef Sum of number of derefined tetrahedra across all chares.
     502                 :            : //! \param[in] sumrefmode Sum of contributions from all chares, encoding
     503                 :            : //!   refinement mode of operation.
     504                 :            : // *****************************************************************************
     505                 :            : {
     506                 :         38 :   auto meshid = tk::cref_find( m_meshid, summeshid );
     507                 :            : 
     508                 :            :   // If at least a single edge on a chare still needs correction, do correction,
     509                 :            :   // otherwise, this mesh refinement step is complete
     510         [ -  + ]:         38 :   if (nextra > 0) {
     511                 :            : 
     512                 :          0 :     ++m_ncit[meshid];
     513                 :          0 :     m_refiner[meshid].comExtra();
     514                 :            : 
     515                 :            :   } else {
     516                 :            : 
     517                 :         38 :     tk::Print print;
     518                 :            : 
     519                 :            :     // decode refmode
     520                 :            :     auto refmode = static_cast< Refiner::RefMode >(
     521                 :         38 :                      sumrefmode / static_cast<std::size_t>(m_nchare[meshid]) );
     522                 :            : 
     523         [ +  - ]:         38 :     if (refmode == Refiner::RefMode::T0REF) {
     524                 :            : 
     525         [ +  + ]:         38 :       if (!g_cfg.get< tag::feedback >()) {
     526                 :         35 :         const auto& initref = g_cfg.get< tag::href_init >();
     527         [ +  - ]:         35 :         print << '\n';
     528 [ +  - ][ +  - ]:        525 :         print.diag( { "meshid", "t0ref", "type", "nref", "nderef", "ncorr" },
         [ +  - ][ +  + ]
         [ +  + ][ -  - ]
                 [ -  - ]
     529                 :            :                     { std::to_string(meshid),
     530                 :         35 :                       std::to_string(m_nt0refit[meshid]),
     531                 :         35 :                       initref[ m_nt0refit[ meshid ] ],
     532                 :            :                       std::to_string(nref),
     533                 :            :                       std::to_string(nderef),
     534                 :         35 :                       std::to_string(m_ncit[meshid]) } );
     535                 :         35 :         ++m_nt0refit[meshid];
     536 [ +  + ][ +  - ]:         35 :         if (m_nt0refit[meshid] == initref.size()) print << '\n';
     537                 :            :       }
     538         [ +  - ]:         38 :       m_progMesh.inc< REFINE >( print );
     539                 :            : 
     540         [ -  - ]:          0 :     } else if (refmode == Refiner::RefMode::DTREF) {
     541                 :            : 
     542 [ -  - ][ -  - ]:          0 :       print.diag( { "meshid", "dtref", "type", "nref", "nderef", "ncorr" },
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
                 [ -  - ]
     543                 :            :                   { std::to_string(meshid),
     544                 :          0 :                     std::to_string(++m_ndtrefit[meshid]),
     545                 :            :                     "error",
     546                 :            :                     std::to_string(nref),
     547                 :            :                     std::to_string(nderef),
     548                 :          0 :                     std::to_string(m_ncit[meshid]) } );
     549                 :            : 
     550 [ -  - ][ -  - ]:          0 :     } else Throw( "RefMode not implemented" );
                 [ -  - ]
     551                 :            : 
     552                 :         38 :     m_ncit[meshid] = 0;
     553         [ +  - ]:         38 :     m_refiner[meshid].perform();
     554                 :            : 
     555                 :            :   }
     556 [ +  - ][ +  - ]:        178 : }
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
     557                 :            : 
     558                 :            : void
     559                 :         38 : Transporter::bndint( tk::real sx, tk::real sy, tk::real sz, tk::real cb,
     560                 :            :                      tk::real summeshid )
     561                 :            : // *****************************************************************************
     562                 :            : // Compute surface integral across the whole problem and perform leak-test
     563                 :            : //! \param[in] sx X component of vector summed
     564                 :            : //! \param[in] sy Y component of vector summed
     565                 :            : //! \param[in] sz Z component of vector summed
     566                 :            : //! \param[in] cb Invoke callback if positive
     567                 :            : //! \param[in] summeshid Mesh id (summed accross all chares)
     568                 :            : //! \details This function aggregates partial surface integrals across the
     569                 :            : //!   boundary faces of the whole problem. After this global sum a
     570                 :            : //!   non-zero vector result indicates a leak, e.g., a hole in the boundary,
     571                 :            : //!   which indicates an error in the boundary face data structures used to
     572                 :            : //!   compute the partial surface integrals.
     573                 :            : // *****************************************************************************
     574                 :            : {
     575         [ +  - ]:         38 :   /*auto meshid =*/tk::cref_find( m_meshid, static_cast<std::size_t>(summeshid) );
     576                 :            : 
     577         [ +  - ]:         38 :   std::stringstream err;
     578         [ +  - ]:         38 :   if (cb < 0.0) {
     579                 :            :     err << "Mesh boundary leaky after mesh refinement step; this is due to a "
     580                 :            :      "problem with updating the side sets used to specify boundary conditions "
     581         [ +  - ]:         38 :      "on faces: ";
     582         [ -  - ]:          0 :   } else if (cb > 0.0) {
     583                 :            :     err << "Mesh boundary leaky during initialization; this is due to "
     584                 :            :     "incorrect or incompletely specified boundary conditions for a given input "
     585         [ -  - ]:          0 :     "mesh: ";
     586                 :            :   }
     587                 :            : 
     588                 :         38 :   auto eps = 1.0e-10;
     589 [ +  - ][ +  - ]:         38 :   if (std::abs(sx) > eps || std::abs(sy) > eps || std::abs(sz) > eps) {
         [ -  + ][ -  + ]
     590         [ -  - ]:          0 :     err << "Integral result must be a zero vector: " << std::setprecision(12) <<
     591 [ -  - ][ -  - ]:          0 :            std::abs(sx) << ", " << std::abs(sy) << ", " << std::abs(sz) <<
         [ -  - ][ -  - ]
                 [ -  - ]
     592 [ -  - ][ -  - ]:          0 :            ", eps = " << eps;
     593 [ -  - ][ -  - ]:          0 :     Throw( err.str() );
                 [ -  - ]
     594                 :            :   }
     595                 :         38 : }
     596                 :            : 
     597                 :            : void
     598                 :        253 : Transporter::refined( std::size_t summeshid,
     599                 :            :                       std::size_t nelem,
     600                 :            :                       std::size_t npoin )
     601                 :            : // *****************************************************************************
     602                 :            : // Reduction target: all chares have refined their mesh
     603                 :            : //! \param[in] summeshid Mesh id (summed accross all Refiner chares)
     604                 :            : //! \param[in] nelem Total number of elements in mesh summed across the
     605                 :            : //!   distributed mesh
     606                 :            : //! \param[in] npoin Total number of mesh points summed across the distributed
     607                 :            : //!   mesh. Note that in parallel this is larger than the number of points in
     608                 :            : //!   the mesh, because the boundary nodes are multi-counted. But we only need
     609                 :            : //!   an equal or larger than npoin for Sorter::setup, so this is okay.
     610                 :            : // *****************************************************************************
     611                 :            : {
     612                 :        253 :   auto meshid = tk::cref_find( m_meshid, summeshid );
     613                 :            : 
     614                 :            :   // Store new number of elements for initially refined mesh
     615                 :        253 :   m_nelem[meshid] = nelem;
     616                 :            : 
     617                 :        253 :   m_sorter[meshid].doneInserting();
     618                 :        253 :   m_sorter[meshid].setup( npoin );
     619                 :        253 : }
     620                 :            : 
     621                 :            : void
     622                 :        253 : Transporter::queriedPart( std::size_t meshid )
     623                 :            : // *****************************************************************************
     624                 :            : // Reduction target: all Partitioner nodes have queried their mesh graphs
     625                 :            : //! \param[in] meshid Mesh id
     626                 :            : // *****************************************************************************
     627                 :            : {
     628                 :        253 :   m_partitioner[meshid].response();
     629                 :        253 : }
     630                 :            : 
     631                 :            : void
     632                 :        253 : Transporter::respondedPart( std::size_t meshid )
     633                 :            : // *****************************************************************************
     634                 :            : // Reduction target: all Partitioner nodes have responded with their mesh graphs
     635                 :            : //! \param[in] meshid Mesh id
     636                 :            : // *****************************************************************************
     637                 :            : {
     638                 :        253 :   m_partitioner[meshid].load();
     639                 :        253 : }
     640                 :            : 
     641                 :            : void
     642                 :        253 : Transporter::queried( std::size_t meshid )
     643                 :            : // *****************************************************************************
     644                 :            : // Reduction target: all Sorter chares have queried their boundary edges
     645                 :            : //! \param[in] meshid Mesh id
     646                 :            : // *****************************************************************************
     647                 :            : {
     648                 :        253 :   m_sorter[meshid].response();
     649                 :        253 : }
     650                 :            : 
     651                 :            : void
     652                 :        253 : Transporter::responded( std::size_t meshid )
     653                 :            : // *****************************************************************************
     654                 :            : // Reduction target: all Sorter chares have responded with their boundary edges
     655                 :            : //! \param[in] meshid Mesh id
     656                 :            : // *****************************************************************************
     657                 :            : {
     658                 :        253 :   m_sorter[meshid].start();
     659                 :        253 : }
     660                 :            : 
     661                 :            : void
     662                 :          0 : Transporter::resized( std::size_t meshid )
     663                 :            : // *****************************************************************************
     664                 :            : // Reduction target: all worker chares have resized their own mesh data after
     665                 :            : //! \param[in] meshid Mesh id
     666                 :            : //! \note Only used for nodal schemes
     667                 :            : // *****************************************************************************
     668                 :            : {
     669                 :          0 :   m_discretization[ meshid ].vol();
     670                 :            : 
     671                 :          0 :   const auto& solver = g_cfg.get< tag::solver >();
     672         [ -  - ]:          0 :   if (solver == "riecg") {
     673                 :          0 :     m_riecg[ meshid ].feop();
     674                 :            :   }
     675         [ -  - ]:          0 :   else if (solver == "laxcg") {
     676                 :          0 :     m_laxcg[ meshid ].feop();
     677                 :            :   }
     678         [ -  - ]:          0 :   else if (solver == "zalcg") {
     679                 :          0 :     m_zalcg[ meshid ].feop();
     680                 :            :   }
     681         [ -  - ]:          0 :   else if (solver == "kozcg") {
     682                 :          0 :     m_kozcg[ meshid ].feop();
     683                 :            :   }
     684         [ -  - ]:          0 :   else if (solver == "chocg") {
     685                 :          0 :     m_chocg[ meshid ].feop();
     686                 :            :   }
     687         [ -  - ]:          0 :   else if (solver == "lohcg") {
     688                 :          0 :     m_lohcg[ meshid ].feop();
     689                 :            :   }
     690                 :            :   else {
     691 [ -  - ][ -  - ]:          0 :     Throw( "Unknown solver: " + solver  );
                 [ -  - ]
     692                 :            :   }
     693                 :          0 : }
     694                 :            : 
     695                 :            : void
     696                 :        253 : Transporter::discinserted( std::size_t meshid )
     697                 :            : // *****************************************************************************
     698                 :            : // Reduction target: all Discretization chares have been inserted
     699                 :            : //! \param[in] meshid Mesh id
     700                 :            : // *****************************************************************************
     701                 :            : {
     702                 :        253 :   m_discretization[ meshid ].doneInserting();
     703                 :        253 : }
     704                 :            : 
     705                 :            : void
     706                 :        265 : Transporter::meshstat( const std::string& header ) const
     707                 :            : // *****************************************************************************
     708                 :            : // Print out mesh statistics
     709                 :            : //! \param[in] header Section header
     710                 :            : // *****************************************************************************
     711                 :            : {
     712                 :        265 :   tk::Print print;
     713                 :            : 
     714         [ +  - ]:        265 :   print.section( header );
     715                 :            : 
     716         [ -  + ]:        265 :   if (m_nelem.size() > 1) {
     717 [ -  - ][ -  - ]:          0 :     print.item( "Number of tetrahedra (per mesh)",tk::parameters(m_nelem) );
                 [ -  - ]
     718 [ -  - ][ -  - ]:          0 :     print.item( "Number of points (per mesh)", tk::parameters(m_npoin) );
                 [ -  - ]
     719 [ -  - ][ -  - ]:          0 :     print.item( "Number of work units (per mesh)", tk::parameters(m_nchare) );
                 [ -  - ]
     720                 :            :   }
     721                 :            : 
     722 [ +  - ][ +  - ]:        265 :   print.item( "Total number of tetrahedra",
     723                 :        265 :               std::accumulate( begin(m_nelem), end(m_nelem), 0UL ) );
     724 [ +  - ][ +  - ]:        265 :   print.item( "Total number of points",
     725                 :        265 :               std::accumulate( begin(m_npoin), end(m_npoin), 0UL ) );
     726 [ +  - ][ +  - ]:        265 :   print.item( "Total number of work units",
     727                 :        265 :               std::accumulate( begin(m_nchare), end(m_nchare), 0 ) );
     728                 :        265 : }
     729                 :            : 
     730                 :            : void
     731                 :        253 : Transporter::disccreated( std::size_t summeshid, std::size_t npoin )
     732                 :            : // *****************************************************************************
     733                 :            : // Reduction target: all Discretization constructors have been called
     734                 :            : //! \param[in] summeshid Mesh id (summed accross all chares)
     735                 :            : //! \param[in] npoin Total number of mesh points (summed across all chares)
     736                 :            : //!  Note that as opposed to npoin in refined(), this npoin is not
     737                 :            : //!  multi-counted, and thus should be correct in parallel.
     738                 :            : // *****************************************************************************
     739                 :            : {
     740                 :        253 :   auto meshid = tk::cref_find( m_meshid, summeshid );
     741                 :            : 
     742                 :            :   // Update number of mesh points for mesh, since it may have been refined
     743         [ +  + ]:        253 :   if (g_cfg.get< tag::href_t0 >()) m_npoin[meshid] = npoin;
     744                 :            : 
     745         [ +  - ]:        253 :   if (++m_ndisc == m_nelem.size()) { // all Disc arrays have been created
     746                 :        253 :     m_ndisc = 0;
     747                 :        253 :     tk::Print print;
     748         [ +  - ]:        253 :     m_progMesh.end( print );
     749         [ +  + ]:        253 :     if (g_cfg.get< tag::href_t0 >()) {
     750 [ +  - ][ +  - ]:         12 :       meshstat( "Mesh initially refined" );
     751                 :            :     }
     752                 :            :   }
     753                 :            : 
     754                 :        253 :   m_refiner[ meshid ].sendProxy();
     755                 :        253 :   m_discretization[ meshid ].vol();
     756                 :            : 
     757         [ +  - ]:        253 :   m_discretization[0][0].npoin(
     758         [ +  - ]:        253 :     std::accumulate( begin(m_npoin), end(m_npoin), 0UL ) );
     759                 :        253 : }
     760                 :            : 
     761                 :            : void
     762                 :        253 : Transporter::workinserted( std::size_t meshid )
     763                 :            : // *****************************************************************************
     764                 :            : // Reduction target: all worker (derived discretization) chares have been
     765                 :            : // inserted
     766                 :            : //! \param[in] meshid Mesh id
     767                 :            : // *****************************************************************************
     768                 :            : {
     769                 :        253 :   const auto& solver = g_cfg.get< tag::solver >();
     770         [ +  + ]:        253 :   if (solver == "riecg") {
     771                 :         73 :     m_riecg[ meshid ].doneInserting();
     772                 :            :   }
     773         [ +  + ]:        180 :   else if (solver == "laxcg") {
     774                 :         23 :     m_laxcg[ meshid ].doneInserting();
     775                 :            :   }
     776         [ +  + ]:        157 :   else if (solver == "zalcg") {
     777                 :         36 :     m_zalcg[ meshid ].doneInserting();
     778                 :            :   }
     779         [ +  + ]:        121 :   else if (solver == "kozcg") {
     780                 :         44 :     m_kozcg[ meshid ].doneInserting();
     781                 :            :   }
     782         [ +  + ]:         77 :   else if (solver == "chocg") {
     783                 :         44 :     m_chocg[ meshid ].doneInserting();
     784                 :         44 :     m_cgpre[ meshid ].doneInserting();
     785                 :         44 :     m_cgmom[ meshid ].doneInserting();
     786                 :            :   }
     787         [ +  - ]:         33 :   else if (solver == "lohcg") {
     788                 :         33 :     m_lohcg[ meshid ].doneInserting();
     789                 :         33 :     m_cgpre[ meshid ].doneInserting();
     790                 :            :   }
     791                 :            :   else {
     792 [ -  - ][ -  - ]:          0 :     Throw( "Unknown solver: " + solver );
                 [ -  - ]
     793                 :            :   }
     794                 :        253 : }
     795                 :            : 
     796                 :            : void
     797                 :        254 : Transporter::diagHeader()
     798                 :            : // *****************************************************************************
     799                 :            : // Configure and write diagnostics file header
     800                 :            : // *****************************************************************************
     801                 :            : {
     802                 :            :   // Output header for diagnostics output file
     803                 :        254 :   tk::DiagWriter dw( g_cfg.get< tag::diag >(),
     804                 :        254 :                      g_cfg.get< tag::diag_format >(),
     805         [ +  - ]:        508 :                      g_cfg.get< tag::diag_precision >() );
     806                 :            : 
     807                 :        254 :   std::vector< std::string > d;
     808                 :            : 
     809                 :        254 :   const auto& solver = g_cfg.get< tag::solver >();
     810         [ +  + ]:        435 :   if (solver == "riecg" ||
     811         [ +  + ]:        339 :       solver == "laxcg" ||
     812 [ +  + ][ +  + ]:        593 :       solver == "zalcg" ||
                 [ +  + ]
     813                 :        122 :       solver == "kozcg")
     814                 :            :   {
     815                 :            : 
     816                 :            :     // Collect variables names for integral/diagnostics output
     817 [ +  - ][ +  + ]:       1056 :     std::vector< std::string > var{ "r", "ru", "rv", "rw", "rE" };
                 [ -  - ]
     818                 :        352 :     auto ncomp = g_cfg.get< tag::problem_ncomp >();
     819         [ +  + ]:        200 :     for (std::size_t c=5; c<ncomp; ++c)
     820 [ +  - ][ +  - ]:         24 :       var.push_back( "c" + std::to_string(c-5) );
                 [ +  - ]
     821                 :            : 
     822                 :        176 :     auto nv = var.size();
     823                 :            : 
     824                 :            :     // Add 'L2(var)' for all variables
     825 [ +  - ][ +  - ]:       1080 :     for (std::size_t i=0; i<nv; ++i) d.push_back( "L2(" + var[i] + ')' );
         [ +  - ][ +  + ]
     826                 :            : 
     827                 :            :     // Add L2-norm of the residuals
     828 [ +  - ][ +  - ]:       1080 :     for (std::size_t i=0; i<nv; ++i) d.push_back( "L2(d" + var[i] + ')' );
         [ +  - ][ +  + ]
     829                 :            : 
     830                 :            :     // Add total energy
     831 [ +  - ][ +  - ]:        176 :     d.push_back( "mE" );
     832                 :            : 
     833                 :            :     // Augment diagnostics variables by error norms (if computed)
     834 [ +  - ][ +  + ]:        176 :     if (problems::SOL()) {
     835 [ +  - ][ +  - ]:         62 :       d.push_back( "L2(err:r)" );
     836 [ +  - ][ +  - ]:         62 :       d.push_back( "L2(err:u)" );
     837 [ +  - ][ +  - ]:         62 :       d.push_back( "L2(err:v)" );
     838 [ +  - ][ +  - ]:         62 :       d.push_back( "L2(err:w)" );
     839 [ +  - ][ +  - ]:         62 :       d.push_back( "L2(err:e)" );
     840 [ +  - ][ +  - ]:         84 :       for (std::size_t i=5; i<nv; ++i) d.push_back( "L2(err:" + var[i] + ')' );
         [ +  - ][ +  + ]
     841 [ +  - ][ +  - ]:         62 :       d.push_back( "L1(err:r)" );
     842 [ +  - ][ +  - ]:         62 :       d.push_back( "L1(err:u)" );
     843 [ +  - ][ +  - ]:         62 :       d.push_back( "L1(err:v)" );
     844 [ +  - ][ +  - ]:         62 :       d.push_back( "L1(err:w)" );
     845 [ +  - ][ +  - ]:         62 :       d.push_back( "L1(err:e)" );
     846 [ +  - ][ +  - ]:         84 :       for (std::size_t i=5; i<nv; ++i) d.push_back( "L1(err:" + var[i] + ')' );
         [ +  - ][ +  + ]
     847                 :            :     }
     848                 :            : 
     849                 :        176 :   }
     850         [ +  + ]:         78 :   else if (solver == "chocg") {
     851                 :            : 
     852                 :            :     // query function to evaluate analytic solution (if defined)
     853         [ +  - ]:         45 :     auto pressure_sol = problems::PRESSURE_SOL();
     854                 :            : 
     855                 :            :     // Collect variables names for integral/diagnostics output
     856 [ +  - ][ +  + ]:         90 :     std::vector< std::string > var{ "p" };
                 [ -  - ]
     857         [ +  + ]:         45 :     if (!pressure_sol) {
     858 [ +  - ][ +  - ]:         38 :       var.push_back( "u" );
     859 [ +  - ][ +  - ]:         38 :       var.push_back( "v" );
     860 [ +  - ][ +  - ]:         38 :       var.push_back( "w" );
     861                 :            :     }
     862                 :            : 
     863                 :         45 :     auto ncomp = g_cfg.get< tag::problem_ncomp >();
     864         [ +  + ]:         49 :     for (std::size_t c=3; c<ncomp; ++c) {
     865 [ +  - ][ +  - ]:          4 :       var.push_back( "c" + std::to_string(c-3) );
                 [ +  - ]
     866                 :            :     }
     867                 :            : 
     868                 :         45 :     auto nv = var.size();
     869                 :            : 
     870                 :            :     // Add 'L2(var)' for all variables
     871 [ +  - ][ +  - ]:        208 :     for (std::size_t i=0; i<nv; ++i) d.push_back( "L2(" + var[i] + ')' );
         [ +  - ][ +  + ]
     872                 :            : 
     873                 :            :     // Add L2-norm of the residuals
     874 [ +  - ][ +  - ]:        208 :     for (std::size_t i=0; i<nv; ++i) d.push_back( "L2(d" + var[i] + ')' );
         [ +  - ][ +  + ]
     875                 :            : 
     876                 :            :     // Augment diagnostics variables by error norms of pressure (if computed)
     877         [ +  + ]:         45 :     if (pressure_sol) {
     878 [ +  - ][ +  - ]:          7 :       d.push_back( "L2(err:p)" );
     879 [ +  - ][ +  - ]:          7 :       d.push_back( "L1(err:p)" );
     880                 :            :     }
     881                 :            :     // Augment diagnostics variables by error norms of adv/diff (if computed)
     882 [ +  - ][ +  + ]:         38 :     else if (problems::SOL()) {
     883 [ +  - ][ +  - ]:          4 :       d.push_back( "L2(err:u)" );
     884 [ +  - ][ +  - ]:          4 :       d.push_back( "L2(err:v)" );
     885 [ +  - ][ +  - ]:          4 :       d.push_back( "L2(err:w)" );
     886 [ +  - ][ +  - ]:          8 :       for (std::size_t i=4; i<nv; ++i) d.push_back( "L2(err:" + var[i] + ')' );
         [ +  - ][ +  + ]
     887 [ +  - ][ +  - ]:          4 :       d.push_back( "L1(err:u)" );
     888 [ +  - ][ +  - ]:          4 :       d.push_back( "L1(err:v)" );
     889 [ +  - ][ +  - ]:          4 :       d.push_back( "L1(err:w)" );
     890 [ +  - ][ +  - ]:          8 :       for (std::size_t i=4; i<nv; ++i) d.push_back( "L1(err:" + var[i] + ')' );
         [ +  - ][ +  + ]
     891                 :            :     }
     892                 :            : 
     893                 :         45 :   }
     894         [ +  - ]:         33 :   else if (solver == "lohcg") {
     895                 :            : 
     896                 :            :     // Collect variables names for integral/diagnostics output
     897 [ +  - ][ +  + ]:         66 :     std::vector< std::string > var{ "p" };
                 [ -  - ]
     898 [ +  - ][ +  - ]:         33 :     var.push_back( "u" );
     899 [ +  - ][ +  - ]:         33 :     var.push_back( "v" );
     900 [ +  - ][ +  - ]:         33 :     var.push_back( "w" );
     901                 :            : 
     902                 :         33 :     auto ncomp = g_cfg.get< tag::problem_ncomp >();
     903         [ +  + ]:         37 :     for (std::size_t c=4; c<ncomp; ++c) {
     904 [ +  - ][ +  - ]:          4 :       var.push_back( "c" + std::to_string(c-4) );
                 [ +  - ]
     905                 :            :     }
     906                 :            : 
     907                 :         33 :     auto nv = var.size();
     908                 :            : 
     909                 :            :     // Add 'L2(var)' for all variables
     910 [ +  - ][ +  - ]:        169 :     for (std::size_t i=0; i<nv; ++i) d.push_back( "L2(" + var[i] + ')' );
         [ +  - ][ +  + ]
     911                 :            : 
     912                 :            :     // Add L2-norm of the residuals
     913 [ +  - ][ +  - ]:        169 :     for (std::size_t i=0; i<nv; ++i) d.push_back( "L2(d" + var[i] + ')' );
         [ +  - ][ +  + ]
     914                 :            : 
     915                 :            :     // Augment diagnostics variables by error norms of adv/diff (if computed)
     916 [ +  - ][ +  + ]:         33 :     if (problems::SOL()) {
     917 [ +  - ][ +  - ]:          4 :       d.push_back( "L2(err:u)" );
     918 [ +  - ][ +  - ]:          4 :       d.push_back( "L2(err:v)" );
     919 [ +  - ][ +  - ]:          4 :       d.push_back( "L2(err:w)" );
     920 [ +  - ][ +  - ]:          8 :       for (std::size_t i=4; i<nv; ++i) d.push_back( "L2(err:" + var[i] + ')' );
         [ +  - ][ +  + ]
     921 [ +  - ][ +  - ]:          4 :       d.push_back( "L1(err:u)" );
     922 [ +  - ][ +  - ]:          4 :       d.push_back( "L1(err:v)" );
     923 [ +  - ][ +  - ]:          4 :       d.push_back( "L1(err:w)" );
     924 [ +  - ][ +  - ]:          8 :       for (std::size_t i=4; i<nv; ++i) d.push_back( "L1(err:" + var[i] + ')' );
         [ +  - ][ +  + ]
     925                 :            :     }
     926                 :            : 
     927                 :         33 :   }
     928                 :            :   else {
     929 [ -  - ][ -  - ]:          0 :     Throw( "Unknown solver: " + solver );
                 [ -  - ]
     930                 :            :   }
     931                 :            : 
     932                 :            :   // Write diagnostics header
     933         [ +  - ]:        254 :   dw.header( d );
     934 [ +  - ][ +  - ]:        762 : }
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ -  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
                 [ -  - ]
     935                 :            : 
     936                 :            : void
     937                 :        254 : Transporter::integralsHeader()
     938                 :            : // *****************************************************************************
     939                 :            : // Configure and write integrals file header
     940                 :            : // *****************************************************************************
     941                 :            : {
     942                 :        254 :   const auto& sidesets_integral = g_cfg.get< tag::integout >();
     943                 :            : 
     944         [ +  + ]:        254 :   if (sidesets_integral.empty()) return;
     945                 :            : 
     946         [ +  - ]:          8 :   auto filename = g_cfg.get< tag::output >() + ".int";
     947                 :            :   tk::DiagWriter dw( filename,
     948                 :          8 :                      g_cfg.get< tag::integout_format >(),
     949         [ +  - ]:          8 :                      g_cfg.get< tag::integout_precision >() );
     950                 :            : 
     951                 :            :   // Collect variables names for integral output
     952                 :          8 :   std::vector< std::string > var;
     953                 :          8 :   const auto& reqv = g_cfg.get< tag::integout_integrals >();
     954         [ +  - ]:          8 :   std::unordered_set< std::string > req( begin(reqv), end(reqv) );
     955         [ +  + ]:         21 :   for (auto s : sidesets_integral) {
     956 [ +  - ][ +  - ]:         13 :     if (req.count( "mass_flow_rate" )) {
                 [ +  + ]
     957 [ +  - ][ +  - ]:         10 :       var.push_back( "mass_flow_rate:" + std::to_string(s) );
     958                 :            :     }
     959 [ +  - ][ +  - ]:         13 :     if (req.count( "force" )) {
                 [ +  + ]
     960                 :          3 :       auto si = std::to_string( s );
     961 [ +  - ][ +  - ]:          3 :       var.push_back( "force_x:" + si );
     962 [ +  - ][ +  - ]:          3 :       var.push_back( "force_y:" + si );
     963 [ +  - ][ +  - ]:          3 :       var.push_back( "force_z:" + si );
     964                 :          3 :     }
     965                 :            :   }
     966                 :            : 
     967                 :            :   // Write integrals header
     968         [ +  - ]:          8 :   dw.header( var );
     969                 :          8 : }
     970                 :            : 
     971                 :            : void
     972                 :        253 : Transporter::totalvol( tk::real v, tk::real initial, tk::real summeshid )
     973                 :            : // *****************************************************************************
     974                 :            : // Reduction target summing total mesh volume across all workers
     975                 :            : //! \param[in] v Mesh volume summed across the distributed mesh
     976                 :            : //! \param[in] initial Sum of contributions from all chares. If larger than
     977                 :            : //!    zero, we are during setup, if zero, during time stepping.
     978                 :            : //! \param[in] summeshid Mesh id (summed accross the distributed mesh)
     979                 :            : // *****************************************************************************
     980                 :            : {
     981         [ +  - ]:        253 :   auto meshid = tk::cref_find( m_meshid, static_cast<std::size_t>(summeshid) );
     982                 :            : 
     983                 :        253 :   m_meshvol[meshid] = v;
     984                 :            : 
     985         [ +  - ]:        253 :   if (initial > 0.0) {   // during initialization
     986                 :            : 
     987                 :        253 :     m_discretization[ meshid ].stat( v );
     988                 :            : 
     989                 :            :   } else {               // during AMR
     990                 :            : 
     991                 :          0 :     const auto& solver = g_cfg.get< tag::solver >();
     992         [ -  - ]:          0 :     if (solver == "riecg") {
     993                 :          0 :       m_riecg[ meshid ].resize_complete();
     994                 :            :     }
     995         [ -  - ]:          0 :     else if (solver == "laxcg") {
     996                 :          0 :       m_laxcg[ meshid ].resize_complete();
     997                 :            :     }
     998         [ -  - ]:          0 :     else if (solver == "zalcg") {
     999                 :          0 :       m_zalcg[ meshid ].resize_complete();
    1000                 :            :     }
    1001         [ -  - ]:          0 :     else if (solver == "kozcg") {
    1002                 :          0 :       m_kozcg[ meshid ].resize_complete();
    1003                 :            :     }
    1004         [ -  - ]:          0 :     else if (solver == "chocg") {
    1005                 :          0 :       m_chocg[ meshid ].resize_complete();
    1006                 :            :     }
    1007         [ -  - ]:          0 :     else if (solver == "lohcg") {
    1008                 :          0 :       m_lohcg[ meshid ].resize_complete();
    1009                 :            :     }
    1010                 :            :     else {
    1011 [ -  - ][ -  - ]:          0 :       Throw( "Unknown solver: " + solver );
                 [ -  - ]
    1012                 :            :     }
    1013                 :            : 
    1014                 :            :   }
    1015                 :        253 : }
    1016                 :            : 
    1017                 :            : void
    1018                 :        253 : Transporter::minstat( tk::real d0, tk::real d1, tk::real d2, tk::real d3,
    1019                 :            :                       tk::real d4, tk::real d5, tk::real rmeshid )
    1020                 :            : // *****************************************************************************
    1021                 :            : // Reduction target yielding minimum mesh statistcs across all workers
    1022                 :            : //! \param[in] d0 Minimum mesh statistics collected over all chares
    1023                 :            : //! \param[in] d1 Minimum mesh statistics collected over all chares
    1024                 :            : //! \param[in] d2 Minimum mesh statistics collected over all chares
    1025                 :            : //! \param[in] d3 Minimum mesh statistics collected over all chares
    1026                 :            : //! \param[in] d4 Minimum mesh statistics collected over all chares
    1027                 :            : //! \param[in] d5 Minimum mesh statistics collected over all chares
    1028                 :            : //! \param[in] rmeshid Mesh id as a real
    1029                 :            : // *****************************************************************************
    1030                 :            : {
    1031                 :        253 :   auto meshid = static_cast<std::size_t>(rmeshid);
    1032                 :            : 
    1033                 :        253 :   m_minstat[meshid][0] = d0;  // minimum edge length
    1034                 :        253 :   m_minstat[meshid][1] = d1;  // minimum cell volume cubic root
    1035                 :        253 :   m_minstat[meshid][2] = d2;  // minimum number of elements on chare
    1036                 :        253 :   m_minstat[meshid][3] = d3;  // minimum number of points on chare
    1037                 :        253 :   m_minstat[meshid][4] = d4;  // minimum number of edges on chare
    1038                 :        253 :   m_minstat[meshid][5] = d5;  // minimum number of comm/total points on chare
    1039                 :            : 
    1040                 :        253 :   minstat_complete(meshid);
    1041                 :        253 : }
    1042                 :            : 
    1043                 :            : void
    1044                 :        253 : Transporter::maxstat( tk::real d0, tk::real d1, tk::real d2, tk::real d3,
    1045                 :            :                       tk::real d4, tk::real d5, tk::real rmeshid )
    1046                 :            : // *****************************************************************************
    1047                 :            : // Reduction target yielding the maximum mesh statistics across all workers
    1048                 :            : //! \param[in] d0 Maximum mesh statistics collected over all chares
    1049                 :            : //! \param[in] d1 Maximum mesh statistics collected over all chares
    1050                 :            : //! \param[in] d2 Maximum mesh statistics collected over all chares
    1051                 :            : //! \param[in] d3 Maximum mesh statistics collected over all chares
    1052                 :            : //! \param[in] d4 Minimum mesh statistics collected over all chares
    1053                 :            : //! \param[in] d5 Minimum mesh statistics collected over all chares
    1054                 :            : //! \param[in] rmeshid Mesh id as a real
    1055                 :            : // *****************************************************************************
    1056                 :            : {
    1057                 :        253 :   auto meshid = static_cast<std::size_t>(rmeshid);
    1058                 :            : 
    1059                 :        253 :   m_maxstat[meshid][0] = d0;  // maximum edge length
    1060                 :        253 :   m_maxstat[meshid][1] = d1;  // maximum cell volume cubic root
    1061                 :        253 :   m_maxstat[meshid][2] = d2;  // maximum number of elements on chare
    1062                 :        253 :   m_maxstat[meshid][3] = d3;  // maximum number of points on chare
    1063                 :        253 :   m_maxstat[meshid][4] = d4;  // maximum number of edges on chare
    1064                 :        253 :   m_maxstat[meshid][5] = d5;  // maximum number of comm/total points on chare
    1065                 :            : 
    1066                 :        253 :   maxstat_complete(meshid);
    1067                 :        253 : }
    1068                 :            : 
    1069                 :            : void
    1070                 :        253 : Transporter::sumstat( tk::real d0, tk::real d1, tk::real d2, tk::real d3,
    1071                 :            :                       tk::real d4, tk::real d5, tk::real d6, tk::real d7,
    1072                 :            :                       tk::real d8, tk::real summeshid )
    1073                 :            : // *****************************************************************************
    1074                 :            : // Reduction target yielding the sum mesh statistics across all workers
    1075                 :            : //! \param[in] d0 Sum mesh statistics collected over all chares
    1076                 :            : //! \param[in] d1 Sum mesh statistics collected over all chares
    1077                 :            : //! \param[in] d2 Sum mesh statistics collected over all chares
    1078                 :            : //! \param[in] d3 Sum mesh statistics collected over all chares
    1079                 :            : //! \param[in] d4 Sum mesh statistics collected over all chares
    1080                 :            : //! \param[in] d5 Sum mesh statistics collected over all chares
    1081                 :            : //! \param[in] d6 Sum mesh statistics collected over all chares
    1082                 :            : //! \param[in] d7 Sum mesh statistics collected over all chares
    1083                 :            : //! \param[in] d8 Sum mesh statistics collected over all chares
    1084                 :            : //! \param[in] summeshid Mesh id (summed accross the distributed mesh)
    1085                 :            : // *****************************************************************************
    1086                 :            : {
    1087         [ +  - ]:        253 :   auto meshid = tk::cref_find( m_meshid, static_cast<std::size_t>(summeshid) );
    1088                 :            : 
    1089                 :        253 :   m_avgstat[meshid][0] = d1 / d0;  // avg edge length
    1090                 :        253 :   m_avgstat[meshid][1] = d3 / d2;  // avg cell volume cubic root
    1091                 :        253 :   m_avgstat[meshid][2] = d5 / d4;  // avg number of elements per chare
    1092                 :        253 :   m_avgstat[meshid][3] = d6 / d4;  // avg number of points per chare
    1093                 :        253 :   m_avgstat[meshid][4] = d7 / d4;  // avg number of edges per chare
    1094                 :        253 :   m_avgstat[meshid][5] = d8 / d4;  // avg number of comm/total points per chare
    1095                 :            : 
    1096                 :        253 :   sumstat_complete(meshid);
    1097                 :        253 : }
    1098                 :            : 
    1099                 :            : void
    1100                 :        253 : Transporter::pdfstat( CkReductionMsg* msg )
    1101                 :            : // *****************************************************************************
    1102                 :            : // Reduction target yielding PDF of mesh statistics across all workers
    1103                 :            : //! \param[in] msg Serialized PDF
    1104                 :            : // *****************************************************************************
    1105                 :            : {
    1106                 :            :   std::size_t meshid;
    1107                 :        253 :   std::vector< tk::UniPDF > pdf;
    1108                 :            : 
    1109                 :            :   // Deserialize final PDF
    1110                 :        253 :   PUP::fromMem creator( msg->getData() );
    1111                 :            :   // cppcheck-suppress uninitvar
    1112         [ +  - ]:        253 :   creator | meshid;
    1113         [ +  - ]:        253 :   creator | pdf;
    1114 [ +  - ][ +  - ]:        253 :   delete msg;
    1115                 :            : 
    1116                 :            :   // cppcheck-suppress uninitvar
    1117         [ +  - ]:        253 :   auto id = std::to_string(meshid);
    1118                 :            : 
    1119                 :            :   // Create new PDF file (overwrite if exists)
    1120 [ +  - ][ +  - ]:        253 :   tk::PDFWriter pdfe( "mesh_edge_pdf." + id + ".txt" );
         [ +  - ][ +  - ]
    1121                 :            :   // Output edgelength PDF
    1122                 :            :   // cppcheck-suppress containerOutOfBounds
    1123         [ +  - ]:        253 :   pdfe.writeTxt( pdf[0],
    1124                 :        506 :                  tk::ctr::PDFInfo{ {"PDF"}, {}, {"edgelength"}, 0, 0.0 } );
    1125                 :            : 
    1126                 :            :   // Create new PDF file (overwrite if exists)
    1127 [ +  - ][ +  - ]:        253 :   tk::PDFWriter pdfv( "mesh_vol_pdf." + id + ".txt" );
         [ +  - ][ +  - ]
    1128                 :            :   // Output cell volume cubic root PDF
    1129                 :            :   // cppcheck-suppress containerOutOfBounds
    1130         [ +  - ]:        253 :   pdfv.writeTxt( pdf[1],
    1131                 :        506 :                  tk::ctr::PDFInfo{ {"PDF"}, {}, {"V^{1/3}"}, 0, 0.0 } );
    1132                 :            : 
    1133                 :            :   // Create new PDF file (overwrite if exists)
    1134 [ +  - ][ +  - ]:        253 :   tk::PDFWriter pdfn( "mesh_ntet_pdf." + id + ".txt" );
         [ +  - ][ +  - ]
    1135                 :            :   // Output number of cells PDF
    1136                 :            :   // cppcheck-suppress containerOutOfBounds
    1137         [ +  - ]:        253 :   pdfn.writeTxt( pdf[2],
    1138                 :        506 :                  tk::ctr::PDFInfo{ {"PDF"}, {}, {"ntets"}, 0, 0.0 } );
    1139                 :            : 
    1140                 :        253 :   pdfstat_complete(meshid);
    1141 [ +  - ][ +  - ]:       2277 : }
         [ +  - ][ +  + ]
         [ +  - ][ +  - ]
         [ +  - ][ +  + ]
         [ +  - ][ +  - ]
         [ +  - ][ +  + ]
         [ +  - ][ -  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
                 [ -  - ]
    1142                 :            : 
    1143                 :            : void
    1144                 :        253 : Transporter::stat()
    1145                 :            : // *****************************************************************************
    1146                 :            : // Echo diagnostics on mesh statistics
    1147                 :            : // *****************************************************************************
    1148                 :            : {
    1149                 :        253 :   tk::Print print;
    1150                 :            : 
    1151         [ +  - ]:        253 :   if (++m_nstat == m_nelem.size()) {     // stats from all meshes have arrived
    1152                 :        253 :     m_nstat = 0;
    1153                 :            : 
    1154         [ +  - ]:        253 :     auto& t = tk::ref_find( m_timer, TimerTag::MESH_PART );
    1155         [ +  - ]:        253 :     print << '\n';
    1156 [ +  - ][ +  - ]:        253 :     print << "Mesh partitioning time: " + std::to_string(t.second) + " sec\n";
         [ +  - ][ +  - ]
    1157         [ +  - ]:        253 :     t = tk::ref_find( m_timer, TimerTag::MESH_DIST );
    1158 [ +  - ][ +  - ]:        253 :     print << "Mesh distribution time: " + std::to_string(t.second) + " sec\n";
         [ +  - ][ +  - ]
    1159                 :            : 
    1160                 :            : 
    1161         [ +  + ]:        506 :     for (std::size_t i=0; i<m_nelem.size(); ++i) {
    1162         [ -  + ]:        253 :       if (m_nelem.size() > 1) {
    1163 [ -  - ][ -  - ]:          0 :         print.section("Mesh " + std::to_string(i) + " distribution statistics");
         [ -  - ][ -  - ]
    1164                 :            :       } else {
    1165 [ +  - ][ +  - ]:        253 :         print.section( "Mesh distribution statistics" );
    1166                 :            :       }
    1167                 :            :       print <<
    1168         [ +  - ]:        506 :         "min/max/avg edgelength = " +
    1169 [ +  - ][ +  - ]:       1012 :         std::to_string( m_minstat[i][0] ) + " / " +
                 [ +  - ]
    1170 [ +  - ][ +  - ]:       1012 :         std::to_string( m_maxstat[i][0] ) + " / " +
                 [ +  - ]
    1171 [ +  - ][ +  - ]:       1012 :         std::to_string( m_avgstat[i][0] ) + "\n" +
                 [ +  - ]
    1172         [ +  - ]:        506 :         "min/max/avg V^(1/3) = " +
    1173 [ +  - ][ +  - ]:       1012 :         std::to_string( m_minstat[i][1] ) + " / " +
                 [ +  - ]
    1174 [ +  - ][ +  - ]:       1012 :         std::to_string( m_maxstat[i][1] ) + " / " +
                 [ +  - ]
    1175 [ +  - ][ +  - ]:       1012 :         std::to_string( m_avgstat[i][1] ) + "\n" +
                 [ +  - ]
    1176         [ +  - ]:        506 :         "min/max/avg nelem = " +
    1177 [ +  - ][ +  - ]:       1012 :         std::to_string( static_cast<std::size_t>(m_minstat[i][2]) ) + " / " +
                 [ +  - ]
    1178 [ +  - ][ +  - ]:       1012 :         std::to_string( static_cast<std::size_t>(m_maxstat[i][2]) ) + " / " +
                 [ +  - ]
    1179 [ +  - ][ +  - ]:       1012 :         std::to_string( static_cast<std::size_t>(m_avgstat[i][2]) ) + "\n" +
                 [ +  - ]
    1180         [ +  - ]:        506 :         "min/max/avg npoin = " +
    1181 [ +  - ][ +  - ]:       1012 :         std::to_string( static_cast<std::size_t>(m_minstat[i][3]) ) + " / " +
                 [ +  - ]
    1182 [ +  - ][ +  - ]:       1012 :         std::to_string( static_cast<std::size_t>(m_maxstat[i][3]) ) + " / " +
                 [ +  - ]
    1183 [ +  - ][ +  - ]:       1012 :         std::to_string( static_cast<std::size_t>(m_avgstat[i][3]) ) + "\n" +
                 [ +  - ]
    1184         [ +  - ]:        506 :         "min/max/avg nedge = " +
    1185 [ +  - ][ +  - ]:       1012 :         std::to_string( static_cast<std::size_t>(m_minstat[i][4]) ) + " / " +
                 [ +  - ]
    1186 [ +  - ][ +  - ]:       1012 :         std::to_string( static_cast<std::size_t>(m_maxstat[i][4]) ) + " / " +
                 [ +  - ]
    1187 [ +  - ][ +  - ]:       1012 :         std::to_string( static_cast<std::size_t>(m_avgstat[i][4]) ) + '\n' +
                 [ +  - ]
    1188         [ +  - ]:        506 :         "min/max/avg ncompoin/npoin = " +
    1189 [ +  - ][ +  - ]:       1012 :         std::to_string( m_minstat[i][5] ) + " / " +
                 [ +  - ]
    1190 [ +  - ][ +  - ]:       1012 :         std::to_string( m_maxstat[i][5] ) + " / " +
                 [ +  - ]
    1191 [ +  - ][ +  - ]:        759 :         std::to_string( m_avgstat[i][5] ) + '\n';
                 [ +  - ]
    1192                 :            :     }
    1193                 :            : 
    1194                 :            :     // Print out time integration header to screen
    1195         [ +  - ]:        253 :     inthead( print );
    1196                 :            : 
    1197 [ +  - ][ +  - ]:        253 :     m_progWork.start( print, "Preparing workers", {{ m_nchare[0] }} );
    1198                 :            : 
    1199                 :            :     // Create "derived-class" workers
    1200 [ +  - ][ +  + ]:        506 :     for (std::size_t i=0; i<m_nelem.size(); ++i) m_sorter[i].createWorkers();
    1201                 :            :   }
    1202                 :        253 : }
    1203                 :            : 
    1204                 :            : void
    1205                 :        253 : Transporter::boxvol( tk::real v, tk::real summeshid )
    1206                 :            : // *****************************************************************************
    1207                 :            : // Reduction target computing total volume of IC box(es)
    1208                 :            : //! \param[in] v Total volume within user-specified IC box(es)
    1209                 :            : //! \param[in] summeshid Mesh id as a real (summed accross the distributed mesh)
    1210                 :            : // *****************************************************************************
    1211                 :            : {
    1212         [ +  - ]:        253 :   auto meshid = tk::cref_find( m_meshid, static_cast<std::size_t>(summeshid) );
    1213 [ +  + ][ +  - ]:        253 :   if (v > 0.0) tk::Print() << "IC-box-volume sum: " + std::to_string(v) << '\n';
         [ +  - ][ +  - ]
                 [ +  - ]
    1214                 :            : 
    1215                 :        253 :   const auto& solver = g_cfg.get< tag::solver >();
    1216         [ +  + ]:        253 :   if (solver == "riecg") {
    1217                 :         73 :     m_riecg[ meshid ].setup( v );
    1218                 :            :   }
    1219         [ +  + ]:        180 :   else if (solver == "laxcg") {
    1220                 :         23 :     m_laxcg[ meshid ].setup( v );
    1221                 :            :   }
    1222         [ +  + ]:        157 :   else if (solver == "zalcg") {
    1223                 :         36 :     m_zalcg[ meshid ].setup( v );
    1224                 :            :   }
    1225         [ +  + ]:        121 :   else if (solver == "kozcg") {
    1226                 :         44 :     m_kozcg[ meshid ].setup( v );
    1227                 :            :   }
    1228         [ +  + ]:         77 :   else if (solver == "chocg") {
    1229                 :         44 :     m_chocg[ meshid ].setup( v );
    1230                 :            :   }
    1231         [ +  - ]:         33 :   else if (solver == "lohcg") {
    1232                 :         33 :     m_lohcg[ meshid ].setup( v );
    1233                 :            :   }
    1234                 :            :   else {
    1235 [ -  - ][ -  - ]:          0 :     Throw( "Unknown solver: " + solver );
                 [ -  - ]
    1236                 :            :   }
    1237                 :            : 
    1238                 :            :   // Turn on automatic load balancing
    1239         [ +  - ]:        253 :   if (++m_ncom == m_nelem.size()) { // all worker arrays have finished
    1240                 :        253 :     m_ncom = 0;
    1241                 :        253 :     tk::Print print;
    1242         [ +  - ]:        253 :     m_progWork.end( print );
    1243         [ +  - ]:        253 :     tk::CProxy_LBSwitch::ckNew();
    1244                 :            :   }
    1245                 :        253 : }
    1246                 :            : 
    1247                 :            : void
    1248                 :        265 : Transporter::inthead( const tk::Print& print )
    1249                 :            : // *****************************************************************************
    1250                 :            : // Print out time integration header to screen
    1251                 :            : //! \param[in] print Pretty printer object to use for printing
    1252                 :            : // *****************************************************************************
    1253                 :            : {
    1254                 :        265 :   const auto dea = g_cfg.get< tag::deactivate >();
    1255         [ +  - ]:        265 :   const auto solver = g_cfg.get< tag::solver >();
    1256         [ +  + ]:        265 :   const auto pre = solver == "chocg" ? 1 : 0;
    1257                 :        265 :   const auto theta = g_cfg.get< tag::theta >();
    1258                 :        265 :   const auto eps = std::numeric_limits< tk::real >::epsilon();
    1259 [ +  + ][ +  + ]:        265 :   const auto mom = solver == "chocg" and theta > eps ? 1 : 0;
    1260                 :            : 
    1261 [ +  - ][ +  - ]:        265 :   print.section( "Time integration" );
    1262                 :            :   print <<
    1263                 :            :   "Legend: it - iteration count\n"
    1264                 :            :   "         t - physics time\n"
    1265                 :            :   "        dt - physics time step size\n"
    1266                 :            :   "       ETE - estimated wall-clock time elapsed (h:m:s)\n"
    1267                 :            :   "       ETA - estimated wall-clock time for accomplishment (h:m:s)\n"
    1268                 :            :   "       EGT - estimated grind wall-clock time (1e-6sec/timestep)\n"
    1269                 :            :   "       EGP - estimated grind performance: wall-clock time "
    1270                 :            :                 "(1e-6sec/DOF/timestep)\n"
    1271                 :            :   "       flg - status flags, legend:\n"
    1272                 :            :   "             f - field (volume and surface) output\n"
    1273                 :            :   "             i - integral output\n"
    1274                 :            :   "             d - diagnostics output\n"
    1275                 :            :   "             t - physics time history output\n"
    1276                 :            :   "             h - h-refinement\n"
    1277                 :            :   "             l - load balancing\n"
    1278                 :        265 :   "             c - checkpoint\n" << (dea ?
    1279                 :        265 :   "             e:x/y - x of y work units deactivated\n" : "") << (pre ?
    1280                 :        265 :   "             p:it - pressure linear solve iterations\n" : "") << (mom ?
    1281                 :            :   "             m:it - momentum/transport linear solve iterations\n" : "") <<
    1282                 :            :   "\n      it             t            dt        ETE        ETA        EGT"
    1283                 :            :   "           EGP  flg\n"
    1284                 :            :   "-----------------------------------------------------------------------"
    1285 [ +  - ][ +  + ]:       1060 :   "-----------------\n";
         [ +  - ][ +  + ]
         [ +  - ][ +  + ]
         [ +  - ][ +  - ]
    1286                 :        265 : }
    1287                 :            : 
    1288                 :            : void
    1289                 :       3201 : Transporter::rhodiagnostics( CkReductionMsg* msg )
    1290                 :            : // *****************************************************************************
    1291                 :            : //  Reduction target collecting diagnostics from density-based solvers
    1292                 :            : //! \param[in] msg Serialized diagnostics vector aggregated across all PEs
    1293                 :            : // *****************************************************************************
    1294                 :            : {
    1295                 :            :   using namespace diagnostics;
    1296                 :            : 
    1297                 :            :   std::size_t meshid;
    1298                 :            :   std::size_t ncomp;
    1299                 :       3201 :   std::vector< std::vector< tk::real > > d;
    1300                 :            : 
    1301                 :            :   // Deserialize diagnostics vector
    1302                 :       3201 :   PUP::fromMem creator( msg->getData() );
    1303                 :            :   // cppcheck-suppress uninitvar
    1304         [ +  - ]:       3201 :   creator | meshid;
    1305         [ +  - ]:       3201 :   creator | ncomp;
    1306         [ +  - ]:       3201 :   creator | d;
    1307 [ +  - ][ +  - ]:       3201 :   delete msg;
    1308                 :            : 
    1309                 :            :   // cppcheck-suppress uninitvar
    1310                 :            :   // cppcheck-suppress unreadVariable
    1311         [ +  - ]:       3201 :   auto id = std::to_string(meshid);
    1312                 :            : 
    1313 [ -  + ][ -  - ]:       3201 :   Assert( ncomp > 0, "Number of scalar components must be positive");
         [ -  - ][ -  - ]
    1314 [ -  + ][ -  - ]:       3201 :   Assert( d.size() == NUMDIAG, "Diagnostics vector size mismatch" );
         [ -  - ][ -  - ]
    1315                 :            : 
    1316                 :            :   // cppcheck-suppress unsignedLessThanZero
    1317         [ +  + ]:      28809 :   for (std::size_t i=0; i<d.size(); ++i) {
    1318 [ -  + ][ -  - ]:      25608 :      Assert( d[i].size() == ncomp, "Size mismatch at final stage of "
         [ -  - ][ -  - ]
    1319                 :            :              "diagnostics aggregation for mesh " + id );
    1320                 :            :   }
    1321                 :            : 
    1322                 :            :   // Allocate storage for those diagnostics that are always computed
    1323         [ +  - ]:       3201 :   std::vector< tk::real > diag( ncomp, 0.0 );
    1324                 :            : 
    1325                 :            :   // Finish computing the L2 norm of conserved variables
    1326         [ +  + ]:      19480 :   for (std::size_t i=0; i<d[L2SOL].size(); ++i) {
    1327                 :            :     // cppcheck-suppress uninitvar
    1328                 :      16279 :     diag[i] = sqrt( d[L2SOL][i] / m_meshvol[meshid] );
    1329                 :            :   }
    1330                 :            :  
    1331                 :            :   // Finish computing the L2 norm of the residuals
    1332         [ +  - ]:       3201 :   std::vector< tk::real > l2res( d[L2RES].size(), 0.0 );
    1333         [ +  + ]:      19480 :   for (std::size_t i=0; i<d[L2RES].size(); ++i) {
    1334                 :            :     // cppcheck-suppress uninitvar
    1335                 :      16279 :     l2res[i] = std::sqrt( d[L2RES][i] / m_meshvol[meshid] );
    1336         [ +  - ]:      16279 :     diag.push_back( l2res[i] );
    1337                 :            :   }
    1338                 :            : 
    1339                 :            :   // Append total energy
    1340         [ +  - ]:       3201 :   diag.push_back( d[TOTALEN][0] );
    1341                 :            : 
    1342                 :            :   // Finish computing norms of the numerical - analytical solution
    1343 [ +  - ][ +  + ]:       3201 :   if (problems::SOL()) {
    1344         [ +  + ]:      10200 :     for (std::size_t i=0; i<d[L2ERR].size(); ++i) {
    1345                 :            :       // cppcheck-suppress uninitvar
    1346         [ +  - ]:       8544 :       diag.push_back( std::sqrt( d[L2ERR][i] / m_meshvol[meshid] ) );
    1347                 :            :     }
    1348         [ +  + ]:      10200 :     for (std::size_t i=0; i<d[L1ERR].size(); ++i) {
    1349                 :            :       // cppcheck-suppress uninitvar
    1350         [ +  - ]:       8544 :       diag.push_back( d[L1ERR][i] / m_meshvol[meshid] );
    1351                 :            :     }
    1352                 :            :   }
    1353                 :            :  
    1354                 :            :   // Append diagnostics file at selected times
    1355         [ +  - ]:       3201 :   auto filename = g_cfg.get< tag::diag >();
    1356 [ -  + ][ -  - ]:       3201 :   if (m_nelem.size() > 1) filename += '.' + id;
                 [ -  - ]
    1357                 :            :   tk::DiagWriter dw( filename,
    1358                 :       3201 :                      g_cfg.get< tag::diag_format >(),
    1359                 :       3201 :                      g_cfg.get< tag::diag_precision >(),
    1360         [ +  - ]:       3201 :                      std::ios_base::app );
    1361         [ +  - ]:       3201 :   dw.write( static_cast<uint64_t>(d[ITER][0]), d[TIME][0], d[DT][0], diag );
    1362                 :            : 
    1363                 :       3201 :   const auto& solver = g_cfg.get< tag::solver >();
    1364         [ +  + ]:       3201 :   if (solver == "riecg") {
    1365                 :            :     // cppcheck-suppress uninitvar
    1366         [ +  - ]:       1613 :     m_riecg[ meshid ].evalres( l2res );
    1367                 :            :   }
    1368         [ +  + ]:       1588 :   else if (solver == "laxcg") {
    1369                 :            :     // cppcheck-suppress uninitvar
    1370         [ +  - ]:        260 :     m_laxcg[ meshid ].evalres( l2res );
    1371                 :            :   }
    1372         [ +  + ]:       1328 :   else if (solver == "zalcg") {
    1373                 :            :     // cppcheck-suppress uninitvar
    1374         [ +  - ]:        480 :     m_zalcg[ meshid ].evalres( l2res );
    1375                 :            :   }
    1376         [ +  - ]:        848 :   else if (solver == "kozcg") {
    1377                 :            :     // cppcheck-suppress uninitvar
    1378         [ +  - ]:        848 :     m_kozcg[ meshid ].evalres( l2res );
    1379                 :            :   }
    1380                 :            :   else {
    1381 [ -  - ][ -  - ]:          0 :     Throw( "Unknown solver: " + solver );
                 [ -  - ]
    1382                 :            :   }
    1383                 :       3201 : }
    1384                 :            : 
    1385                 :            : void
    1386                 :        517 : Transporter::prediagnostics( CkReductionMsg* msg )
    1387                 :            : // *****************************************************************************
    1388                 :            : //  Reduction target collecting diagnostics from pressure-based solvers
    1389                 :            : //! \param[in] msg Serialized diagnostics vector aggregated across all PEs
    1390                 :            : // *****************************************************************************
    1391                 :            : {
    1392                 :            :   using namespace diagnostics;
    1393                 :            : 
    1394                 :            :   std::size_t meshid;
    1395                 :            :   std::size_t ncomp;
    1396                 :        517 :   std::vector< std::vector< tk::real > > d;
    1397                 :            : 
    1398                 :            :   // Deserialize diagnostics vector
    1399                 :        517 :   PUP::fromMem creator( msg->getData() );
    1400                 :            :   // cppcheck-suppress uninitvar
    1401         [ +  - ]:        517 :   creator | meshid;
    1402         [ +  - ]:        517 :   creator | ncomp;
    1403         [ +  - ]:        517 :   creator | d;
    1404 [ +  - ][ +  - ]:        517 :   delete msg;
    1405                 :            : 
    1406                 :            :   // cppcheck-suppress uninitvar
    1407                 :            :   // cppcheck-suppress unreadVariable
    1408         [ +  - ]:        517 :   auto id = std::to_string(meshid);
    1409                 :            : 
    1410 [ -  + ][ -  - ]:        517 :   Assert( ncomp > 0, "Number of scalar components must be positive");
         [ -  - ][ -  - ]
    1411 [ -  + ][ -  - ]:        517 :   Assert( d.size() == NUMDIAG, "Diagnostics vector size mismatch" );
         [ -  - ][ -  - ]
    1412                 :            : 
    1413                 :            :   // cppcheck-suppress unsignedLessThanZero
    1414         [ +  + ]:       4653 :   for (std::size_t i=0; i<d.size(); ++i) {
    1415 [ -  + ][ -  - ]:       4136 :      Assert( d[i].size() == ncomp, "Size mismatch at final stage of "
         [ -  - ][ -  - ]
    1416                 :            :              "diagnostics aggregation for mesh " + id );
    1417                 :            :   }
    1418                 :            : 
    1419                 :            :   // Allocate storage for those diagnostics that are always computed
    1420         [ +  - ]:        517 :   std::vector< tk::real > diag( ncomp, 0.0 );
    1421                 :            : 
    1422                 :            :   // Finish computing the L2 norm of conserved variables
    1423         [ +  + ]:       2644 :   for (std::size_t i=0; i<d[L2SOL].size(); ++i) {
    1424                 :            :     // cppcheck-suppress uninitvar
    1425                 :       2127 :     diag[i] = sqrt( d[L2SOL][i] / m_meshvol[meshid] );
    1426                 :            :   }
    1427                 :            : 
    1428                 :            :   // Finish computing the L2 norm of the residuals
    1429         [ +  - ]:        517 :   std::vector< tk::real > l2res( d[L2RES].size(), 0.0 );
    1430         [ +  + ]:       2644 :   for (std::size_t i=0; i<d[L2RES].size(); ++i) {
    1431                 :            :     // cppcheck-suppress uninitvar
    1432                 :       2127 :     l2res[i] = std::sqrt( d[L2RES][i] / m_meshvol[meshid] );
    1433         [ +  - ]:       2127 :     diag.push_back( l2res[i] );
    1434                 :            :   }
    1435                 :            : 
    1436                 :            :   // Finish computing norms of the numerical - analytical pressure solution
    1437 [ +  - ][ +  + ]:        517 :   if (problems::PRESSURE_SOL()) {
    1438         [ +  - ]:          7 :     diag.push_back( std::sqrt( d[L2ERR][0] / m_meshvol[meshid] ) );
    1439         [ +  - ]:          7 :     diag.push_back( d[L1ERR][0] / m_meshvol[meshid] );
    1440                 :            :   }
    1441                 :            : 
    1442                 :            :   // Finish computing norms of the numerical - analytical adv/diff solution
    1443 [ +  - ][ +  + ]:        517 :   if (problems::SOL()) {
    1444         [ +  + ]:        407 :     for (std::size_t i=1; i<d[L2ERR].size(); ++i) {
    1445                 :            :       // cppcheck-suppress uninitvar
    1446         [ +  - ]:        320 :       diag.push_back( std::sqrt( d[L2ERR][i] / m_meshvol[meshid] ) );
    1447                 :            :     }
    1448         [ +  + ]:        407 :     for (std::size_t i=1; i<d[L1ERR].size(); ++i) {
    1449                 :            :       // cppcheck-suppress uninitvar
    1450         [ +  - ]:        320 :       diag.push_back( d[L1ERR][i] / m_meshvol[meshid] );
    1451                 :            :     }
    1452                 :            :   }
    1453                 :            : 
    1454                 :            :   // Append diagnostics file at selected times
    1455         [ +  - ]:        517 :   auto filename = g_cfg.get< tag::diag >();
    1456 [ -  + ][ -  - ]:        517 :   if (m_nelem.size() > 1) filename += '.' + id;
                 [ -  - ]
    1457                 :            :   tk::DiagWriter dw( filename,
    1458                 :        517 :                      g_cfg.get< tag::diag_format >(),
    1459                 :        517 :                      g_cfg.get< tag::diag_precision >(),
    1460         [ +  - ]:        517 :                      std::ios_base::app );
    1461         [ +  - ]:        517 :   dw.write( static_cast<uint64_t>(d[ITER][0]), d[TIME][0], d[DT][0], diag );
    1462                 :            : 
    1463                 :        517 :   const auto& solver = g_cfg.get< tag::solver >();
    1464         [ +  - ]:        517 :   if (solver == "chocg") {
    1465                 :            :     // cppcheck-suppress uninitvar
    1466         [ +  - ]:        517 :     m_chocg[ meshid ].evalres( l2res );
    1467                 :            :   }
    1468                 :            :   else {
    1469 [ -  - ][ -  - ]:          0 :     Throw( "Unknown solver: " + solver );
                 [ -  - ]
    1470                 :            :   }
    1471                 :        517 : }
    1472                 :            : 
    1473                 :            : void
    1474                 :        424 : Transporter::acdiagnostics( CkReductionMsg* msg )
    1475                 :            : // *****************************************************************************
    1476                 :            : //  Reduction target collecting diagnostics from artificial compressibility
    1477                 :            : //  solvers
    1478                 :            : //! \param[in] msg Serialized diagnostics vector aggregated across all PEs
    1479                 :            : // *****************************************************************************
    1480                 :            : {
    1481                 :            :   using namespace diagnostics;
    1482                 :            : 
    1483                 :            :   std::size_t meshid;
    1484                 :            :   std::size_t ncomp;
    1485                 :        424 :   std::vector< std::vector< tk::real > > d;
    1486                 :            : 
    1487                 :            :   // Deserialize diagnostics vector
    1488                 :        424 :   PUP::fromMem creator( msg->getData() );
    1489                 :            :   // cppcheck-suppress uninitvar
    1490         [ +  - ]:        424 :   creator | meshid;
    1491         [ +  - ]:        424 :   creator | ncomp;
    1492         [ +  - ]:        424 :   creator | d;
    1493 [ +  - ][ +  - ]:        424 :   delete msg;
    1494                 :            : 
    1495                 :            :   // cppcheck-suppress uninitvar
    1496                 :            :   // cppcheck-suppress unreadVariable
    1497         [ +  - ]:        424 :   auto id = std::to_string(meshid);
    1498                 :            : 
    1499 [ -  + ][ -  - ]:        424 :   Assert( ncomp > 0, "Number of scalar components must be positive");
         [ -  - ][ -  - ]
    1500 [ -  + ][ -  - ]:        424 :   Assert( d.size() == NUMDIAG, "Diagnostics vector size mismatch" );
         [ -  - ][ -  - ]
    1501                 :            : 
    1502                 :            :   // cppcheck-suppress unsignedLessThanZero
    1503         [ +  + ]:       3816 :   for (std::size_t i=0; i<d.size(); ++i) {
    1504 [ -  + ][ -  - ]:       3392 :      Assert( d[i].size() == ncomp, "Size mismatch at final stage of "
         [ -  - ][ -  - ]
    1505                 :            :              "diagnostics aggregation for mesh " + id );
    1506                 :            :   }
    1507                 :            : 
    1508                 :            :   // Allocate storage for those diagnostics that are always computed
    1509         [ +  - ]:        424 :   std::vector< tk::real > diag( ncomp, 0.0 );
    1510                 :            : 
    1511                 :            :   // Finish computing the L2 norm of conserved variables
    1512         [ +  + ]:       2200 :   for (std::size_t i=0; i<d[L2SOL].size(); ++i) {
    1513                 :            :     // cppcheck-suppress uninitvar
    1514                 :       1776 :     diag[i] = sqrt( d[L2SOL][i] / m_meshvol[meshid] );
    1515                 :            :   }
    1516                 :            : 
    1517                 :            :   // Finish computing the L2 norm of the residuals
    1518         [ +  - ]:        424 :   std::vector< tk::real > l2res( d[L2RES].size(), 0.0 );
    1519         [ +  + ]:       2200 :   for (std::size_t i=0; i<d[L2RES].size(); ++i) {
    1520                 :            :     // cppcheck-suppress uninitvar
    1521                 :       1776 :     l2res[i] = std::sqrt( d[L2RES][i] / m_meshvol[meshid] );
    1522         [ +  - ]:       1776 :     diag.push_back( l2res[i] );
    1523                 :            :   }
    1524                 :            : 
    1525                 :            :   // Finish computing norms of the numerical - analytical adv/diff solution
    1526 [ +  - ][ +  + ]:        424 :   if (problems::SOL()) {
    1527         [ +  + ]:        400 :     for (std::size_t i=1; i<d[L2ERR].size(); ++i) {
    1528                 :            :       // cppcheck-suppress uninitvar
    1529         [ +  - ]:        320 :       diag.push_back( std::sqrt( d[L2ERR][i] / m_meshvol[meshid] ) );
    1530                 :            :     }
    1531         [ +  + ]:        400 :     for (std::size_t i=1; i<d[L1ERR].size(); ++i) {
    1532                 :            :       // cppcheck-suppress uninitvar
    1533         [ +  - ]:        320 :       diag.push_back( d[L1ERR][i] / m_meshvol[meshid] );
    1534                 :            :     }
    1535                 :            :   }
    1536                 :            : 
    1537                 :            :   // Append diagnostics file at selected times
    1538         [ +  - ]:        424 :   auto filename = g_cfg.get< tag::diag >();
    1539 [ -  + ][ -  - ]:        424 :   if (m_nelem.size() > 1) filename += '.' + id;
                 [ -  - ]
    1540                 :            :   tk::DiagWriter dw( filename,
    1541                 :        424 :                      g_cfg.get< tag::diag_format >(),
    1542                 :        424 :                      g_cfg.get< tag::diag_precision >(),
    1543         [ +  - ]:        424 :                      std::ios_base::app );
    1544         [ +  - ]:        424 :   dw.write( static_cast<uint64_t>(d[ITER][0]), d[TIME][0], d[DT][0], diag );
    1545                 :            : 
    1546                 :        424 :   const auto& solver = g_cfg.get< tag::solver >();
    1547         [ +  - ]:        424 :   if (solver == "lohcg") {
    1548                 :            :     // cppcheck-suppress uninitvar
    1549         [ +  - ]:        424 :     m_lohcg[ meshid ].evalres( l2res );
    1550                 :            :   }
    1551                 :            :   else {
    1552 [ -  - ][ -  - ]:          0 :     Throw( "Unknown solver: " + solver );
                 [ -  - ]
    1553                 :            :   }
    1554                 :        424 : }
    1555                 :            : 
    1556                 :            : void
    1557                 :         98 : Transporter::integrals( CkReductionMsg* msg )
    1558                 :            : // *****************************************************************************
    1559                 :            : // Reduction target optionally collecting integrals
    1560                 :            : //! \param[in] msg Serialized integrals aggregated across all PEs
    1561                 :            : // *****************************************************************************
    1562                 :            : {
    1563                 :            :   using namespace integrals;
    1564                 :            : 
    1565                 :            :   // cppcheck-suppress unassignedVariable
    1566                 :            :   std::size_t meshid;
    1567                 :         98 :   std::vector< std::map< int, tk::real > > d;
    1568                 :            : 
    1569                 :            :   // Deserialize integrals vector
    1570                 :         98 :   PUP::fromMem creator( msg->getData() );
    1571                 :            :   // cppcheck-suppress uninitvar
    1572         [ +  - ]:         98 :   creator | meshid;
    1573         [ +  - ]:         98 :   creator | d;
    1574 [ +  - ][ +  - ]:         98 :   delete msg;
    1575                 :            : 
    1576                 :         98 :   const auto& sidesets_integral = g_cfg.get< tag::integout >();
    1577                 :            :   // cppcheck-suppress
    1578         [ +  - ]:         98 :   if (not sidesets_integral.empty()) {
    1579                 :            : 
    1580 [ -  + ][ -  - ]:         98 :     Assert( d.size() == NUMINT, "Integrals vector size mismatch" );
         [ -  - ][ -  - ]
    1581                 :            : 
    1582                 :            :     // Collect integrals for output
    1583                 :         98 :     std::vector< tk::real > ints;
    1584 [ +  - ][ +  + ]:        246 :     for (const auto& [s,m] : d[MASS_FLOW_RATE]) ints.push_back( m );
    1585 [ +  - ][ +  + ]:        122 :     for (const auto& [s,m] : d[FORCE_X]) ints.push_back( m );
    1586 [ +  - ][ +  + ]:        122 :     for (const auto& [s,m] : d[FORCE_Y]) ints.push_back( m );
    1587 [ +  - ][ +  + ]:        122 :     for (const auto& [s,m] : d[FORCE_Z]) ints.push_back( m );
    1588                 :            : 
    1589                 :            :     // Append integrals file at selected times
    1590         [ +  - ]:         98 :     auto filename = g_cfg.get< tag::output >() + ".int";
    1591                 :            :     tk::DiagWriter dw( filename,
    1592                 :         98 :                        g_cfg.get< tag::integout_format >(),
    1593                 :         98 :                        g_cfg.get< tag::integout_precision >(),
    1594         [ +  - ]:         98 :                        std::ios_base::app );
    1595                 :            :     // cppcheck-suppress containerOutOfBounds
    1596 [ +  - ][ +  - ]:         98 :     dw.write( static_cast<uint64_t>(tk::cref_find( d[ITER], 0 )),
    1597                 :            :               // cppcheck-suppress containerOutOfBounds
    1598         [ +  - ]:         98 :               tk::cref_find( d[TIME], 0 ),
    1599                 :            :               // cppcheck-suppress containerOutOfBounds
    1600         [ +  - ]:         98 :               tk::cref_find( d[DT], 0 ),
    1601                 :            :               ints );
    1602                 :         98 :   }
    1603                 :            : 
    1604                 :         98 :   const auto& solver = g_cfg.get< tag::solver >();
    1605         [ +  + ]:         98 :   if (solver == "riecg") {
    1606                 :            :     // cppcheck-suppress uninitvar
    1607         [ +  - ]:         74 :     m_riecg[ meshid ].step();
    1608                 :            :   }
    1609         [ -  + ]:         24 :   else if (solver == "laxcg") {
    1610                 :            :     // cppcheck-suppress uninitvar
    1611         [ -  - ]:          0 :     m_laxcg[ meshid ].step();
    1612                 :            :   }
    1613         [ -  + ]:         24 :   else if (solver == "zalcg") {
    1614                 :            :     // cppcheck-suppress uninitvar
    1615         [ -  - ]:          0 :     m_zalcg[ meshid ].step();
    1616                 :            :   }
    1617         [ -  + ]:         24 :   else if (solver == "kozcg") {
    1618                 :            :     // cppcheck-suppress uninitvar
    1619         [ -  - ]:          0 :     m_kozcg[ meshid ].step();
    1620                 :            :   }
    1621         [ +  + ]:         24 :   else if (solver == "chocg") {
    1622                 :            :     // cppcheck-suppress uninitvar
    1623         [ +  - ]:         20 :     m_chocg[ meshid ].step();
    1624                 :            :   }
    1625         [ +  - ]:          4 :   else if (solver == "lohcg") {
    1626                 :            :     // cppcheck-suppress uninitvar
    1627         [ +  - ]:          4 :     m_lohcg[ meshid ].step();
    1628                 :            :   }
    1629                 :            :   else
    1630 [ -  - ][ -  - ]:          0 :     Throw( "Unknown solver: " + solver );
                 [ -  - ]
    1631                 :         98 : }
    1632                 :            : 
    1633                 :            : void
    1634                 :        277 : Transporter::resume()
    1635                 :            : // *****************************************************************************
    1636                 :            : // Resume execution from checkpoint/restart files
    1637                 :            : //! \details This is invoked by Charm++ after the checkpoint is done, as well as
    1638                 :            : //!   when the restart (returning from a checkpoint) is complete
    1639                 :            : // *****************************************************************************
    1640                 :            : {
    1641         [ +  + ]:        554 :   if (std::any_of(begin(m_finished), end(m_finished), [](auto f){return !f;})) {
    1642                 :            : 
    1643                 :            :     // If just restarted from a checkpoint, Main( CkMigrateMessage* msg ) has
    1644                 :            :     // increased g_nrestart, but only on PE 0, so broadcast.
    1645                 :            : 
    1646                 :         12 :     const auto& solver = g_cfg.get< tag::solver >();
    1647         [ +  + ]:         12 :     if (solver == "riecg") {
    1648         [ +  + ]:          4 :       for (std::size_t i=0; i<m_nelem.size(); ++i) {
    1649                 :          2 :         m_riecg[i].evalLB( g_nrestart );
    1650                 :            :       }
    1651                 :            :     }
    1652         [ +  + ]:         10 :     else if (solver == "laxcg") {
    1653         [ +  + ]:          4 :       for (std::size_t i=0; i<m_nelem.size(); ++i) {
    1654                 :          2 :         m_laxcg[i].evalLB( g_nrestart );
    1655                 :            :       }
    1656                 :            :     }
    1657         [ +  + ]:          8 :     else if (solver == "zalcg") {
    1658         [ +  + ]:          4 :       for (std::size_t i=0; i<m_nelem.size(); ++i) {
    1659                 :          2 :         m_zalcg[i].evalLB( g_nrestart );
    1660                 :            :       }
    1661                 :            :     }
    1662         [ +  + ]:          6 :     else if ( solver == "kozcg") {
    1663         [ +  + ]:          4 :       for (std::size_t i=0; i<m_nelem.size(); ++i) {
    1664                 :          2 :         m_kozcg[i].evalLB( g_nrestart );
    1665                 :            :       }
    1666                 :            :     }
    1667         [ +  + ]:          4 :     else if ( solver == "chocg") {
    1668         [ +  + ]:          4 :       for (std::size_t i=0; i<m_nelem.size(); ++i) {
    1669                 :          2 :         m_chocg[i].evalLB( g_nrestart );
    1670                 :            :       }
    1671                 :            :     }
    1672         [ +  - ]:          2 :     else if ( solver == "lohcg") {
    1673         [ +  + ]:          4 :       for (std::size_t i=0; i<m_nelem.size(); ++i) {
    1674                 :          2 :         m_lohcg[i].evalLB( g_nrestart );
    1675                 :            :       }
    1676                 :            :     }
    1677                 :            :     else {
    1678 [ -  - ][ -  - ]:          0 :       Throw( "Unknown solver: " + solver );
                 [ -  - ]
    1679                 :            :     }
    1680                 :            : 
    1681                 :            : 
    1682                 :            :   } else {
    1683                 :            : 
    1684                 :        265 :     mainProxy.finalize();
    1685                 :            : 
    1686                 :            :   }
    1687                 :        277 : }
    1688                 :            : 
    1689                 :            : void
    1690                 :        265 : Transporter::checkpoint( std::size_t finished, std::size_t meshid )
    1691                 :            : // *****************************************************************************
    1692                 :            : // Save checkpoint/restart files
    1693                 :            : //! \param[in] finished Nonzero if finished with time stepping
    1694                 :            : //! \param[in] meshid Mesh id
    1695                 :            : // *****************************************************************************
    1696                 :            : {
    1697                 :        265 :   m_finished[meshid] = finished;
    1698                 :            : 
    1699         [ +  - ]:        265 :   if (++m_nchk == m_nelem.size()) { // all worker arrays have checkpointed
    1700                 :        265 :     m_nchk = 0;
    1701         [ +  + ]:        265 :     if (not g_cfg.get< tag::benchmark >()) {
    1702                 :        157 :       const auto& ckptdir = g_cfg.get< tag::checkpoint >();
    1703 [ +  - ][ +  - ]:        157 :       CkCallback res( CkIndex_Transporter::resume(), thisProxy );
    1704         [ +  - ]:        157 :       CkStartCheckpoint( ckptdir.c_str(), res );
    1705                 :            :       //CkStartMemCheckpoint( res );
    1706                 :        157 :     } else {
    1707                 :        108 :       resume();
    1708                 :            :     }
    1709                 :            :   }
    1710                 :        265 : }
    1711                 :            : 
    1712                 :            : void
    1713                 :        265 : Transporter::finish( std::size_t meshid )
    1714                 :            : // *****************************************************************************
    1715                 :            : // Normal finish of time stepping
    1716                 :            : //! \param[in] meshid Mesh id
    1717                 :            : // *****************************************************************************
    1718                 :            : {
    1719                 :        265 :   checkpoint( /* finished = */ 1, meshid );
    1720                 :        265 : }
    1721                 :            : 
    1722                 :            : #include "NoWarning/transporter.def.h"

Generated by: LCOV version 1.16