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 : 248 : Transporter::Transporter() :
53 [ + - ][ - + ]: 744 : m_input{ g_cfg.get< tag::input >() },
[ + + ][ + + ]
[ - - ][ - - ]
54 [ + - ][ + - ]: 248 : m_nchare( m_input.size() ),
55 [ + - ]: 248 : m_ncit( m_nchare.size(), 0 ),
56 : 248 : m_nload( 0 ),
57 : 248 : m_npart( 0 ),
58 : 248 : m_nstat( 0 ),
59 : 248 : m_ndisc( 0 ),
60 : 248 : m_nchk( 0 ),
61 : 248 : m_ncom( 0 ),
62 [ + - ][ - - ]: 248 : m_nt0refit( m_nchare.size(), 0 ),
63 [ + - ][ - - ]: 248 : m_ndtrefit( m_nchare.size(), 0 ),
64 [ + - ][ - - ]: 248 : m_noutrefit( m_nchare.size(), 0 ),
65 [ + - ][ - - ]: 248 : m_noutderefit( m_nchare.size(), 0 ),
66 [ + - ][ + - ]: 248 : m_nelem( m_nchare.size() ),
67 [ + - ][ - - ]: 248 : m_finished( m_nchare.size(), 0 ),
68 [ + - ][ - - ]: 248 : m_meshvol( m_nchare.size() ),
69 [ + - ][ - - ]: 248 : m_minstat( m_nchare.size() ),
70 [ + - ][ - - ]: 248 : m_maxstat( m_nchare.size() ),
71 [ + - ][ + - ]: 248 : m_avgstat( m_nchare.size() ),
[ - - ]
72 [ + - ]: 248 : m_progMesh( g_cfg.get< tag::feedback >(), ProgMeshPrefix, ProgMeshLegend ),
73 [ + - ][ + - ]: 744 : m_progWork( g_cfg.get< tag::feedback >(), ProgWorkPrefix, ProgWorkLegend )
[ + - ]
74 : : // *****************************************************************************
75 : : // Constructor
76 : : // *****************************************************************************
77 : : {
78 : 248 : const auto nstep = g_cfg.get< tag::nstep >();
79 : 248 : const auto t0 = g_cfg.get< tag::t0 >();
80 : 248 : const auto term = g_cfg.get< tag::term >();
81 : 248 : 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 [ + - ][ + - ]: 248 : if ( nstep != 0 && term > t0 && constdt < term-t0 ) {
[ + - ]
90 : :
91 : : // Enable SDAG waits for collecting mesh statistics
92 [ + - ]: 248 : thisProxy.wait4stat();
93 : :
94 : : // Configure and write diagnostics file header
95 [ + - ]: 248 : diagHeader();
96 : :
97 : : // Configure and write integrals file header
98 [ + - ]: 248 : integralsHeader();
99 : :
100 : : // Create mesh partitioner AND boundary condition object group
101 [ + - ]: 248 : createPartitioner();
102 : :
103 [ - - ]: 0 : } else finish(); // stop if no time stepping requested
104 : 247 : }
105 : :
106 : 12 : Transporter::Transporter( CkMigrateMessage* m ) :
107 : : CBase_Transporter( m ),
108 [ + - ]: 12 : m_progMesh( g_cfg.get< tag::feedback >(), ProgMeshPrefix, ProgMeshLegend ),
109 [ + - ][ + - ]: 12 : 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 : : auto print = tk::Print();
116 : : print << "\nXyst> Restarted from checkpoint\n";
117 [ + - ]: 12 : inthead( print );
118 : 12 : }
119 : :
120 : : bool
121 : 247 : 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 : : std::unordered_set< int > usersets;
139 : :
140 : : // Collect side sets at which BCs are set
141 : :
142 [ + + ]: 650 : for (const auto& s : g_cfg.get< tag::bc_dir >()) {
143 [ + - ]: 403 : if (!s.empty()) usersets.insert( s[0] );
144 : : }
145 : :
146 [ + - ][ + + ]: 500 : for (auto s : g_cfg.get< tag::bc_sym >()) usersets.insert( s );
147 : :
148 [ + - ][ + + ]: 257 : for (auto s : g_cfg.get< tag::bc_far >()) usersets.insert( s );
149 : :
150 [ + + ]: 259 : for (const auto& s : g_cfg.get< tag::bc_pre >()) {
151 [ + - ]: 12 : if (!s.empty()) usersets.insert( s[0] );
152 : : }
153 : :
154 [ + + ]: 323 : for (const auto& s : g_cfg.get< tag::pre_bc_dir >()) {
155 [ + - ]: 76 : if (!s.empty()) usersets.insert( s[0] );
156 : : }
157 : :
158 [ + - ][ + + ]: 250 : for (auto s : g_cfg.get< tag::pre_bc_sym >()) usersets.insert( s );
159 : :
160 : : // Add sidesets requested for field output
161 [ + - ][ + + ]: 280 : for (auto s : g_cfg.get< tag::fieldout >()) usersets.insert( s );
162 : : // Add sidesets requested for integral output
163 [ + - ][ + + ]: 260 : 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 : : std::unordered_set< int > sidesets_used;
167 [ + + ]: 996 : for (auto i : usersets) { // for all side sets used in control file
168 [ + - ]: 749 : if (bnd.find(i) != end(bnd)) // used set found among side sets in file
169 : : 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 : 247 : tk::erase_if( bnd, [&]( auto& item ) {
178 [ + + ]: 1376 : return sidesets_used.find( item.first ) == end(sidesets_used);
179 : : });
180 : :
181 : 494 : return not bnd.empty();
182 : : }
183 : :
184 : : void
185 : 248 : Transporter::createPartitioner()
186 : : // *****************************************************************************
187 : : // Create mesh partitioner AND boundary conditions group
188 : : // *****************************************************************************
189 : : {
190 : : // cppcheck-suppress unreadVariable
191 : : auto print = tk::Print();
192 : :
193 : : // Create partitioner callbacks (order important)
194 : : tk::PartitionerCallback cbp {{
195 : 248 : CkCallback( CkReductionTarget(Transporter,queriedPart), thisProxy )
196 : 248 : , CkCallback( CkReductionTarget(Transporter,respondedPart), thisProxy )
197 : 248 : , CkCallback( CkReductionTarget(Transporter,load), thisProxy )
198 : 248 : , CkCallback( CkReductionTarget(Transporter,partitioned), thisProxy )
199 : 248 : , CkCallback( CkReductionTarget(Transporter,distributed), thisProxy )
200 : 248 : , CkCallback( CkReductionTarget(Transporter,refinserted), thisProxy )
201 [ + - ]: 248 : }};
202 : :
203 : : // Create refiner callbacks (order important)
204 : : tk::RefinerCallback cbr {{
205 : 248 : CkCallback( CkReductionTarget(Transporter,queriedRef), thisProxy )
206 : 248 : , CkCallback( CkReductionTarget(Transporter,respondedRef), thisProxy )
207 : 248 : , CkCallback( CkReductionTarget(Transporter,compatibility), thisProxy )
208 : 248 : , CkCallback( CkReductionTarget(Transporter,bndint), thisProxy )
209 : 248 : , CkCallback( CkReductionTarget(Transporter,matched), thisProxy )
210 : 248 : , CkCallback( CkReductionTarget(Transporter,refined), thisProxy )
211 [ + - ]: 248 : }};
212 : :
213 : : // Create sorter callbacks (order important)
214 : : tk::SorterCallback cbs {{
215 : 248 : CkCallback( CkReductionTarget(Transporter,queried), thisProxy )
216 : 248 : , CkCallback( CkReductionTarget(Transporter,responded), thisProxy )
217 : 248 : , CkCallback( CkReductionTarget(Transporter,discinserted), thisProxy )
218 : 248 : , CkCallback( CkReductionTarget(Transporter,workinserted), thisProxy )
219 : 0 : }};
220 : :
221 : : // Start timer measuring preparation of mesh(es) for partitioning
222 [ + - ][ - + ]: 248 : m_timer[ TimerTag::MESH_READ ];
223 : :
224 [ - + ][ - - ]: 248 : ErrChk( !m_input.empty(), "No input mesh" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
225 : :
226 : : // Start preparing mesh(es)
227 [ + - ][ + - ]: 744 : 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 : 248 : std::size_t meshid = 0;
231 [ + + ]: 495 : for (const auto& filename : m_input) {
232 : : // Create mesh reader for reading side sets from file
233 [ + - ]: 248 : tk::ExodusIIMeshReader mr( filename );
234 : :
235 : : // Read out total number of mesh points from mesh file
236 [ + - ][ + - ]: 248 : m_npoin.push_back( mr.npoin() );
[ + - ]
237 : :
238 : : std::map< int, std::vector< std::size_t > > bface;
239 : : std::map< int, std::vector< std::size_t > > faces;
240 : : std::map< int, std::vector< std::size_t > > bnode;
241 : :
242 : : // Read boundary-face connectivity on side sets
243 [ + - ]: 248 : mr.readSidesetFaces( bface, faces );
244 : :
245 : : bool bcs_set = false;
246 : : // Read node lists on side sets
247 [ + - ]: 247 : bnode = mr.readSidesetNodes();
248 : : // Verify boundarty condition (BC) side sets used exist in mesh file
249 [ + - ]: 247 : bcs_set = matchBCs( bnode );
250 [ - + ][ - - ]: 247 : bcs_set = bcs_set || matchBCs( bface );
[ - - ]
251 : :
252 : : // Warn on no BCs
253 : : if (!bcs_set) print << "\n>>> WARNING: No boundary conditions set\n\n";
254 : :
255 : : // Create empty discretization chare array
256 [ + - ][ + - ]: 494 : m_discretization.push_back( CProxy_Discretization::ckNew() );
257 [ + - ]: 247 : CkArrayOptions opt;
258 [ + - ]: 247 : opt.bindTo( m_discretization.back() );
259 : :
260 : : // Create empty discretization scheme chare array (bound to discretization)
261 : : CProxy_RieCG riecg;
262 : : CProxy_LaxCG laxcg;
263 : : CProxy_ZalCG zalcg;
264 : : CProxy_KozCG kozcg;
265 : : CProxy_ChoCG chocg;
266 : : CProxy_LohCG lohcg;
267 : : tk::CProxy_ConjugateGradients cgpre, cgmom;
268 : : const auto& solver = g_cfg.get< tag::solver >();
269 [ + + ]: 247 : if (solver == "riecg") {
270 [ + - ][ + - ]: 146 : m_riecg.push_back( CProxy_RieCG::ckNew(opt) );
[ + - ]
271 : : riecg = m_riecg.back();
272 : : }
273 [ + + ]: 174 : else if (solver == "laxcg") {
274 [ + - ][ + - ]: 46 : m_laxcg.push_back( CProxy_LaxCG::ckNew(opt) );
[ + - ]
275 : : laxcg = m_laxcg.back();
276 : : }
277 [ + + ]: 151 : else if (solver == "zalcg") {
278 [ + - ][ + - ]: 72 : m_zalcg.push_back( CProxy_ZalCG::ckNew(opt) );
[ + - ]
279 : : zalcg = m_zalcg.back();
280 : : }
281 [ + + ]: 115 : else if (solver == "kozcg") {
282 [ + - ][ + - ]: 88 : m_kozcg.push_back( CProxy_KozCG::ckNew(opt) );
[ + - ]
283 : : kozcg = m_kozcg.back();
284 : : }
285 [ + + ]: 71 : else if (solver == "chocg") {
286 [ + - ][ + - ]: 82 : m_chocg.push_back( CProxy_ChoCG::ckNew(opt) );
[ + - ]
287 : : chocg = m_chocg.back();
288 [ + - ][ + - ]: 82 : m_cgpre.push_back( tk::CProxy_ConjugateGradients::ckNew(opt) );
[ + - ]
289 : : cgpre = m_cgpre.back();
290 [ + - ][ + - ]: 82 : m_cgmom.push_back( tk::CProxy_ConjugateGradients::ckNew(opt) );
[ + - ]
291 : : cgmom = m_cgmom.back();
292 : : }
293 [ + - ]: 30 : else if (solver == "lohcg") {
294 [ + - ][ + - ]: 60 : m_lohcg.push_back( CProxy_LohCG::ckNew(opt) );
[ + - ]
295 : : lohcg = m_lohcg.back();
296 [ + - ][ + - ]: 60 : m_cgpre.push_back( tk::CProxy_ConjugateGradients::ckNew(opt) );
[ + - ]
297 : : 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 [ + - ][ + - ]: 494 : m_refiner.push_back( CProxy_Refiner::ckNew(opt) );
305 : : // Create empty mesh sorter Charm++ chare array (bound to discretization)
306 [ + - ][ + - ]: 494 : m_sorter.push_back( CProxy_Sorter::ckNew(opt) );
307 : :
308 : : // Create MeshWriter chare group for mesh
309 [ + - ]: 247 : m_meshwriter.push_back(
310 : : tk::CProxy_MeshWriter::ckNew(
311 [ + - ]: 247 : g_cfg.get< tag::benchmark >(), m_input.size() ) );
312 : :
313 : : // Create mesh partitioner Charm++ chare nodegroup for all meshes
314 [ + - ][ - - ]: 247 : m_partitioner.push_back(
315 : : CProxy_Partitioner::ckNew( meshid, filename, cbp, cbr, cbs,
316 [ + - ]: 247 : thisProxy, m_refiner.back(), m_sorter.back(), m_meshwriter.back(),
317 : : m_discretization.back(), riecg, laxcg, zalcg, kozcg, chocg, lohcg,
318 : : cgpre, cgmom, bface, faces, bnode ) );
319 : :
320 [ - + ]: 247 : ++meshid;
321 : 494 : }
322 : 247 : }
323 : :
324 : : void
325 : 247 : 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 : 247 : meshid /= static_cast< std::size_t >( CkNumNodes() );
334 : : Assert( meshid < m_nelem.size(), "MeshId indexing out" );
335 : 247 : m_nelem[meshid] = nelem;
336 : :
337 : : // Compute load distribution given total work (nelem) and user-specified
338 : : // virtualization
339 : : uint64_t chunksize, remainder;
340 : 247 : m_nchare[meshid] = static_cast<int>(
341 : 247 : tk::linearLoadDistributor(
342 : : g_cfg.get< tag::virt >(),
343 : : 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 : 247 : m_meshid[ static_cast<std::size_t>(m_nchare[meshid])*meshid ] = meshid;
348 : : Assert( meshid < m_nelem.size(), "MeshId indexing out" );
349 : :
350 : : // Partition first mesh
351 [ + - ]: 247 : if (meshid == 0) {
352 : 247 : m_timer[ TimerTag::MESH_PART ]; // start timer measuring mesh partitioning
353 : 247 : m_partitioner[0].partition( m_nchare[0] );
354 : : }
355 : :
356 [ + - ]: 247 : if (++m_nload == m_nelem.size()) { // all meshes have been loaded
357 : 247 : m_nload = 0;
358 : : auto print = tk::Print();
359 : :
360 : : auto& timer = tk::ref_find( m_timer, TimerTag::MESH_READ );
361 : 247 : timer.second = timer.first.dsec();
362 [ + - ]: 247 : print << "Mesh read time: " + std::to_string( timer.second ) + " sec\n";
363 : :
364 : : // Print out mesh partitioning configuration
365 [ + - ]: 247 : print.section( "Partitioning mesh" );
366 [ + - ]: 247 : print.item( "Partitioner", g_cfg.get< tag::part >() );
367 [ + - ]: 247 : print.item( "Virtualization", g_cfg.get< tag::virt >() );
368 : : // Print out initial mesh statistics
369 [ + - ]: 494 : meshstat( "Mesh read from file" );
370 : :
371 : : // Tell meshwriter the total number of chares
372 : 247 : m_meshwriter[meshid].nchare( m_nchare[meshid] );
373 : :
374 : : // Query number of initial mesh refinement steps
375 : : int nref = 0;
376 [ + + ]: 247 : 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 : : int nreord = 0;
382 [ + + ]: 247 : if (g_cfg.get< tag::reorder >()) nreord = m_nchare[0];
383 : :
384 : : print << '\n';
385 : 494 : m_progMesh.start( print, "Preparing mesh", {{ CkNumPes(), CkNumPes(), nref,
386 : : m_nchare[0], m_nchare[0], nreord, nreord }} );
387 : : }
388 : 247 : }
389 : :
390 : : void
391 : 247 : Transporter::partitioned( std::size_t meshid )
392 : : // *****************************************************************************
393 : : // Reduction target: a mesh has been partitioned
394 : : //! \param[in] meshid Mesh id
395 : : // *****************************************************************************
396 : : {
397 [ + - ]: 247 : if (++m_npart == m_nelem.size()) { // all meshes have been partitioned
398 : 247 : m_npart = 0;
399 : : auto& timer = tk::ref_find( m_timer, TimerTag::MESH_PART );
400 : 247 : timer.second = timer.first.dsec();
401 : 247 : 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 : 247 : }
406 : :
407 : : void
408 : 247 : 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 : 247 : m_partitioner[meshid].refine();
416 : : auto& timer = tk::ref_find( m_timer, TimerTag::MESH_DIST );
417 : 247 : timer.second = timer.first.dsec();
418 : 247 : }
419 : :
420 : : void
421 : 247 : 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 [ - + ]: 247 : 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 : : m_refiner[meshid].doneInserting();
447 : :
448 : : }
449 : 247 : }
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 : : 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 : : const auto& initref = g_cfg.get< tag::href_init >();
527 : : print << '\n';
528 [ + - ][ + - ]: 525 : print.diag( { "meshid", "t0ref", "type", "nref", "nderef", "ncorr" },
[ + - ][ + + ]
[ - + ][ + + ]
[ - + ][ + + ]
[ - - ][ - - ]
[ - - ][ - - ]
529 : : { std::to_string(meshid),
530 : : std::to_string(m_nt0refit[meshid]),
531 [ + - ]: 35 : initref[ m_nt0refit[ meshid ] ],
532 : : std::to_string(nref),
533 : : std::to_string(nderef),
534 : : 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 : : std::to_string(++m_ndtrefit[meshid]),
545 : : "error",
546 : : std::to_string(nref),
547 : : std::to_string(nderef),
548 : : 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 : 0 : 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 : 0 : /*auto meshid =*/tk::cref_find( m_meshid, static_cast<std::size_t>(summeshid) );
576 : :
577 : 0 : std::stringstream err;
578 [ - - ]: 0 : 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 [ - - ]: 0 : "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 : : auto eps = 1.0e-10;
589 [ - - ][ - - ]: 0 : if (std::abs(sx) > eps || std::abs(sy) > eps || std::abs(sz) > eps) {
[ - - ]
590 : : err << "Integral result must be a zero vector: " << std::setprecision(12) <<
591 : : std::abs(sx) << ", " << std::abs(sy) << ", " << std::abs(sz) <<
592 : : ", eps = " << eps;
593 [ - - ][ - - ]: 0 : Throw( err.str() );
[ - - ][ - - ]
[ - - ][ - - ]
594 : : }
595 : 0 : }
596 : :
597 : : void
598 : 247 : 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 : 247 : auto meshid = tk::cref_find( m_meshid, summeshid );
613 : :
614 : : // Store new number of elements for initially refined mesh
615 : 247 : m_nelem[meshid] = nelem;
616 : :
617 : : m_sorter[meshid].doneInserting();
618 : 247 : m_sorter[meshid].setup( npoin );
619 : 247 : }
620 : :
621 : : void
622 : 247 : 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 : 247 : m_partitioner[meshid].response();
629 : 247 : }
630 : :
631 : : void
632 : 247 : 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 : 247 : m_partitioner[meshid].load();
639 : 247 : }
640 : :
641 : : void
642 : 247 : 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 : 247 : m_sorter[meshid].response();
649 : 247 : }
650 : :
651 : : void
652 : 247 : 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 : 247 : m_sorter[meshid].start();
659 : 247 : }
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 : : 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 : 247 : 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 : : m_discretization[ meshid ].doneInserting();
703 : 247 : }
704 : :
705 : : void
706 : 259 : Transporter::meshstat( const std::string& header ) const
707 : : // *****************************************************************************
708 : : // Print out mesh statistics
709 : : //! \param[in] header Section header
710 : : // *****************************************************************************
711 : : {
712 : : tk::Print print;
713 : :
714 : 259 : print.section( header );
715 : :
716 [ - + ]: 259 : 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 [ + - ]: 259 : print.item( "Total number of tetrahedra",
723 : 518 : std::accumulate( begin(m_nelem), end(m_nelem), 0UL ) );
724 [ + - ]: 259 : print.item( "Total number of points",
725 : 518 : std::accumulate( begin(m_npoin), end(m_npoin), 0UL ) );
726 [ + - ]: 259 : print.item( "Total number of work units",
727 : 259 : std::accumulate( begin(m_nchare), end(m_nchare), 0 ) );
728 : 259 : }
729 : :
730 : : void
731 [ + + ]: 247 : 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 : 247 : auto meshid = tk::cref_find( m_meshid, summeshid );
741 : :
742 : : // Update number of mesh points for mesh, since it may have been refined
743 [ + + ]: 247 : if (g_cfg.get< tag::href_t0 >()) m_npoin[meshid] = npoin;
744 : :
745 [ + - ]: 247 : if (++m_ndisc == m_nelem.size()) { // all Disc arrays have been created
746 : 247 : m_ndisc = 0;
747 : : tk::Print print;
748 : 247 : m_progMesh.end( print );
749 [ + + ]: 247 : if (g_cfg.get< tag::href_t0 >()) {
750 [ + - ]: 24 : meshstat( "Mesh initially refined" );
751 : : }
752 : : }
753 : :
754 : 247 : m_refiner[ meshid ].sendProxy();
755 : 247 : m_discretization[ meshid ].vol();
756 : :
757 : 247 : m_discretization[0][0].npoin(
758 [ + - ][ - + ]: 247 : std::accumulate( begin(m_npoin), end(m_npoin), 0UL ) );
[ - - ]
759 : 247 : }
760 : :
761 : : void
762 : 247 : 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 : : const auto& solver = g_cfg.get< tag::solver >();
770 [ + + ]: 247 : if (solver == "riecg") {
771 : : m_riecg[ meshid ].doneInserting();
772 : : }
773 [ + + ]: 174 : else if (solver == "laxcg") {
774 : : m_laxcg[ meshid ].doneInserting();
775 : : }
776 [ + + ]: 151 : else if (solver == "zalcg") {
777 : : m_zalcg[ meshid ].doneInserting();
778 : : }
779 [ + + ]: 115 : else if (solver == "kozcg") {
780 : : m_kozcg[ meshid ].doneInserting();
781 : : }
782 [ + + ]: 71 : else if (solver == "chocg") {
783 : : m_chocg[ meshid ].doneInserting();
784 : : m_cgpre[ meshid ].doneInserting();
785 : : m_cgmom[ meshid ].doneInserting();
786 : : }
787 [ + - ]: 30 : else if (solver == "lohcg") {
788 : : m_lohcg[ meshid ].doneInserting();
789 : : m_cgpre[ meshid ].doneInserting();
790 : : }
791 : : else {
792 [ - - ][ - - ]: 0 : Throw( "Unknown solver: " + solver );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
793 : : }
794 : 247 : }
795 : :
796 : : void
797 : 248 : Transporter::diagHeader()
798 : : // *****************************************************************************
799 : : // Configure and write diagnostics file header
800 : : // *****************************************************************************
801 : : {
802 : : // Output header for diagnostics output file
803 : : tk::DiagWriter dw( g_cfg.get< tag::diag >(),
804 : : g_cfg.get< tag::diag_format >(),
805 : 248 : g_cfg.get< tag::diag_precision >() );
806 : :
807 : : std::vector< std::string > d;
808 : :
809 : : const auto& solver = g_cfg.get< tag::solver >();
810 [ + + ]: 175 : if (solver == "riecg" ||
811 [ + + ]: 152 : solver == "laxcg" ||
812 [ + + ][ + + ]: 364 : solver == "zalcg" ||
813 : : 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 : 176 : auto ncomp = g_cfg.get< tag::problem_ncomp >();
819 [ + + ]: 200 : for (std::size_t c=5; c<ncomp; ++c)
820 [ + - ]: 48 : var.push_back( "c" + std::to_string(c-5) );
821 : :
822 : : 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 [ + - ][ + + ]: 352 : 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 [ + + ]: 72 : else if (solver == "chocg") {
851 : :
852 : : // query function to evaluate analytic solution (if defined)
853 [ + - ]: 42 : auto pressure_sol = problems::PRESSURE_SOL();
854 : :
855 : : // Collect variables names for integral/diagnostics output
856 [ - + ][ + + ]: 84 : std::vector< std::string > var{ "p" };
[ - + ][ - - ]
[ - - ]
857 [ + + ]: 42 : if (!pressure_sol) {
858 [ + - ]: 35 : var.push_back( "u" );
859 [ + - ]: 35 : var.push_back( "v" );
860 [ + - ]: 70 : var.push_back( "w" );
861 : : }
862 : :
863 : : auto nv = var.size();
864 : :
865 : : // Add 'L2(var)' for all variables
866 [ + - ][ + - ]: 189 : for (std::size_t i=0; i<nv; ++i) d.push_back( "L2(" + var[i] + ')' );
[ + + ]
867 : :
868 : : // Add L2-norm of the residuals
869 [ + - ][ + - ]: 189 : for (std::size_t i=0; i<nv; ++i) d.push_back( "L2(d" + var[i] + ')' );
[ + + ]
870 : :
871 : : // Augment diagnostics variables by error norms (if computed)
872 [ + + ]: 42 : if (pressure_sol) {
873 [ + - ]: 7 : d.push_back( "L2(err:p)" );
874 [ + - ]: 14 : d.push_back( "L1(err:p)" );
875 : : }
876 : :
877 : 42 : }
878 [ + - ]: 30 : else if (solver == "lohcg") {
879 : :
880 : : // Collect variables names for integral/diagnostics output
881 [ - + ][ + + ]: 60 : std::vector< std::string > var{ "p" };
[ - + ][ - - ]
[ - - ]
882 : 60 : var.push_back( "u" );
883 [ + - ]: 30 : var.push_back( "v" );
884 [ + - ]: 60 : var.push_back( "w" );
885 : :
886 : : auto nv = var.size();
887 : :
888 : : // Add 'L2(var)' for all variables
889 [ + - ][ + - ]: 150 : for (std::size_t i=0; i<nv; ++i) d.push_back( "L2(" + var[i] + ')' );
[ + + ]
890 : :
891 : : // Add L2-norm of the residuals
892 [ + - ][ + - ]: 150 : for (std::size_t i=0; i<nv; ++i) d.push_back( "L2(d" + var[i] + ')' );
[ + + ]
893 : :
894 : 30 : }
895 : : else {
896 [ - - ][ - - ]: 0 : Throw( "Unknown solver: " + solver );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
897 : : }
898 : :
899 : : // Write diagnostics header
900 [ + - ]: 248 : dw.header( d );
901 [ + - ][ + - ]: 992 : }
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ][ + - ]
[ + - ][ - - ]
[ - - ][ - - ]
902 : :
903 : : void
904 [ + + ]: 248 : Transporter::integralsHeader()
905 : : // *****************************************************************************
906 : : // Configure and write integrals file header
907 : : // *****************************************************************************
908 : : {
909 : : const auto& sidesets_integral = g_cfg.get< tag::integout >();
910 : :
911 [ + + ]: 248 : if (sidesets_integral.empty()) return;
912 : :
913 : 8 : auto filename = g_cfg.get< tag::output >() + ".int";
914 : : tk::DiagWriter dw( filename,
915 : : g_cfg.get< tag::integout_format >(),
916 [ + - ]: 8 : g_cfg.get< tag::integout_precision >() );
917 : :
918 : : // Collect variables names for integral output
919 : : std::vector< std::string > var;
920 : : const auto& reqv = g_cfg.get< tag::integout_integrals >();
921 : 8 : std::unordered_set< std::string > req( begin(reqv), end(reqv) );
922 [ + + ]: 21 : for (auto s : sidesets_integral) {
923 [ + - ][ + + ]: 26 : if (req.count( "mass_flow_rate" )) {
924 : 20 : var.push_back( "mass_flow_rate:" + std::to_string(s) );
925 : : }
926 [ + - ][ + + ]: 26 : if (req.count( "force" )) {
927 : 3 : auto si = std::to_string( s );
928 [ + - ]: 3 : var.push_back( "force_x:" + si );
929 [ + - ]: 3 : var.push_back( "force_y:" + si );
930 [ + - ][ - + ]: 6 : var.push_back( "force_z:" + si );
[ - - ]
931 : : }
932 : : }
933 : :
934 : : // Write integrals header
935 [ + - ]: 8 : dw.header( var );
936 : 8 : }
937 : :
938 : : void
939 : 247 : Transporter::totalvol( tk::real v, tk::real initial, tk::real summeshid )
940 : : // *****************************************************************************
941 : : // Reduction target summing total mesh volume across all workers
942 : : //! \param[in] v Mesh volume summed across the distributed mesh
943 : : //! \param[in] initial Sum of contributions from all chares. If larger than
944 : : //! zero, we are during setup, if zero, during time stepping.
945 : : //! \param[in] summeshid Mesh id (summed accross the distributed mesh)
946 : : // *****************************************************************************
947 : : {
948 [ + - ]: 247 : auto meshid = tk::cref_find( m_meshid, static_cast<std::size_t>(summeshid) );
949 : :
950 [ + - ]: 247 : m_meshvol[meshid] = v;
951 : :
952 [ + - ]: 247 : if (initial > 0.0) { // during initialization
953 : :
954 : 247 : m_discretization[ meshid ].stat( v );
955 : :
956 : : } else { // during AMR
957 : :
958 : : const auto& solver = g_cfg.get< tag::solver >();
959 [ - - ]: 0 : if (solver == "riecg") {
960 : 0 : m_riecg[ meshid ].resize_complete();
961 : : }
962 [ - - ]: 0 : else if (solver == "laxcg") {
963 : 0 : m_laxcg[ meshid ].resize_complete();
964 : : }
965 [ - - ]: 0 : else if (solver == "zalcg") {
966 : 0 : m_zalcg[ meshid ].resize_complete();
967 : : }
968 [ - - ]: 0 : else if (solver == "kozcg") {
969 : 0 : m_kozcg[ meshid ].resize_complete();
970 : : }
971 [ - - ]: 0 : else if (solver == "chocg") {
972 : 0 : m_chocg[ meshid ].resize_complete();
973 : : }
974 [ - - ]: 0 : else if (solver == "lohcg") {
975 : 0 : m_lohcg[ meshid ].resize_complete();
976 : : }
977 : : else {
978 [ - - ][ - - ]: 0 : Throw( "Unknown solver: " + solver );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
979 : : }
980 : :
981 : : }
982 : 247 : }
983 : :
984 : : void
985 : 247 : Transporter::minstat( tk::real d0, tk::real d1, tk::real d2, tk::real d3,
986 : : tk::real d4, tk::real d5, tk::real rmeshid )
987 : : // *****************************************************************************
988 : : // Reduction target yielding minimum mesh statistcs across all workers
989 : : //! \param[in] d0 Minimum mesh statistics collected over all chares
990 : : //! \param[in] d1 Minimum mesh statistics collected over all chares
991 : : //! \param[in] d2 Minimum mesh statistics collected over all chares
992 : : //! \param[in] d3 Minimum mesh statistics collected over all chares
993 : : //! \param[in] d4 Minimum mesh statistics collected over all chares
994 : : //! \param[in] d5 Minimum mesh statistics collected over all chares
995 : : //! \param[in] rmeshid Mesh id as a real
996 : : // *****************************************************************************
997 : : {
998 : 247 : auto meshid = static_cast<std::size_t>(rmeshid);
999 : :
1000 : 247 : m_minstat[meshid][0] = d0; // minimum edge length
1001 : 247 : m_minstat[meshid][1] = d1; // minimum cell volume cubic root
1002 : 247 : m_minstat[meshid][2] = d2; // minimum number of elements on chare
1003 : 247 : m_minstat[meshid][3] = d3; // minimum number of points on chare
1004 : 247 : m_minstat[meshid][4] = d4; // minimum number of edges on chare
1005 : 247 : m_minstat[meshid][5] = d5; // minimum number of comm/total points on chare
1006 : :
1007 : 247 : minstat_complete(meshid);
1008 : 247 : }
1009 : :
1010 : : void
1011 : 247 : Transporter::maxstat( tk::real d0, tk::real d1, tk::real d2, tk::real d3,
1012 : : tk::real d4, tk::real d5, tk::real rmeshid )
1013 : : // *****************************************************************************
1014 : : // Reduction target yielding the maximum mesh statistics across all workers
1015 : : //! \param[in] d0 Maximum mesh statistics collected over all chares
1016 : : //! \param[in] d1 Maximum mesh statistics collected over all chares
1017 : : //! \param[in] d2 Maximum mesh statistics collected over all chares
1018 : : //! \param[in] d3 Maximum mesh statistics collected over all chares
1019 : : //! \param[in] d4 Minimum mesh statistics collected over all chares
1020 : : //! \param[in] d5 Minimum mesh statistics collected over all chares
1021 : : //! \param[in] rmeshid Mesh id as a real
1022 : : // *****************************************************************************
1023 : : {
1024 : 247 : auto meshid = static_cast<std::size_t>(rmeshid);
1025 : :
1026 : 247 : m_maxstat[meshid][0] = d0; // maximum edge length
1027 : 247 : m_maxstat[meshid][1] = d1; // maximum cell volume cubic root
1028 : 247 : m_maxstat[meshid][2] = d2; // maximum number of elements on chare
1029 : 247 : m_maxstat[meshid][3] = d3; // maximum number of points on chare
1030 : 247 : m_maxstat[meshid][4] = d4; // maximum number of edges on chare
1031 : 247 : m_maxstat[meshid][5] = d5; // maximum number of comm/total points on chare
1032 : :
1033 : 247 : maxstat_complete(meshid);
1034 : 247 : }
1035 : :
1036 : : void
1037 : 247 : Transporter::sumstat( tk::real d0, tk::real d1, tk::real d2, tk::real d3,
1038 : : tk::real d4, tk::real d5, tk::real d6, tk::real d7,
1039 : : tk::real d8, tk::real summeshid )
1040 : : // *****************************************************************************
1041 : : // Reduction target yielding the sum mesh statistics across all workers
1042 : : //! \param[in] d0 Sum mesh statistics collected over all chares
1043 : : //! \param[in] d1 Sum mesh statistics collected over all chares
1044 : : //! \param[in] d2 Sum mesh statistics collected over all chares
1045 : : //! \param[in] d3 Sum mesh statistics collected over all chares
1046 : : //! \param[in] d4 Sum mesh statistics collected over all chares
1047 : : //! \param[in] d5 Sum mesh statistics collected over all chares
1048 : : //! \param[in] d6 Sum mesh statistics collected over all chares
1049 : : //! \param[in] d7 Sum mesh statistics collected over all chares
1050 : : //! \param[in] d8 Sum mesh statistics collected over all chares
1051 : : //! \param[in] summeshid Mesh id (summed accross the distributed mesh)
1052 : : // *****************************************************************************
1053 : : {
1054 : 247 : auto meshid = tk::cref_find( m_meshid, static_cast<std::size_t>(summeshid) );
1055 : :
1056 : 247 : m_avgstat[meshid][0] = d1 / d0; // avg edge length
1057 : 247 : m_avgstat[meshid][1] = d3 / d2; // avg cell volume cubic root
1058 : 247 : m_avgstat[meshid][2] = d5 / d4; // avg number of elements per chare
1059 : 247 : m_avgstat[meshid][3] = d6 / d4; // avg number of points per chare
1060 : 247 : m_avgstat[meshid][4] = d7 / d4; // avg number of edges per chare
1061 : 247 : m_avgstat[meshid][5] = d8 / d4; // avg number of comm/total points per chare
1062 : :
1063 : 247 : sumstat_complete(meshid);
1064 : 247 : }
1065 : :
1066 : : void
1067 [ + - ]: 247 : Transporter::pdfstat( CkReductionMsg* msg )
1068 : : // *****************************************************************************
1069 : : // Reduction target yielding PDF of mesh statistics across all workers
1070 : : //! \param[in] msg Serialized PDF
1071 : : // *****************************************************************************
1072 : : {
1073 : : std::size_t meshid;
1074 : : std::vector< tk::UniPDF > pdf;
1075 : :
1076 : : // Deserialize final PDF
1077 : : PUP::fromMem creator( msg->getData() );
1078 : : // cppcheck-suppress uninitvar
1079 : : creator | meshid;
1080 : : creator | pdf;
1081 : : delete msg;
1082 : :
1083 : : // cppcheck-suppress uninitvar
1084 [ + - ]: 247 : auto id = std::to_string(meshid);
1085 : :
1086 : : // Create new PDF file (overwrite if exists)
1087 [ + - ][ + - ]: 741 : tk::PDFWriter pdfe( "mesh_edge_pdf." + id + ".txt" );
[ + - ][ - + ]
[ - - ][ - - ]
1088 : : // Output edgelength PDF
1089 : : // cppcheck-suppress containerOutOfBounds
1090 [ + - ]: 247 : pdfe.writeTxt( pdf[0],
1091 [ - + ][ - - ]: 247 : tk::ctr::PDFInfo{ {"PDF"}, {}, {"edgelength"}, 0, 0.0 } );
1092 : :
1093 : : // Create new PDF file (overwrite if exists)
1094 [ + - ][ + - ]: 741 : tk::PDFWriter pdfv( "mesh_vol_pdf." + id + ".txt" );
[ + - ][ - + ]
[ - - ]
1095 : : // Output cell volume cubic root PDF
1096 : : // cppcheck-suppress containerOutOfBounds
1097 [ + - ]: 247 : pdfv.writeTxt( pdf[1],
1098 [ - + ][ - - ]: 247 : tk::ctr::PDFInfo{ {"PDF"}, {}, {"V^{1/3}"}, 0, 0.0 } );
1099 : :
1100 : : // Create new PDF file (overwrite if exists)
1101 [ + - ][ + - ]: 741 : tk::PDFWriter pdfn( "mesh_ntet_pdf." + id + ".txt" );
[ + - ][ - + ]
[ - - ]
1102 : : // Output number of cells PDF
1103 : : // cppcheck-suppress containerOutOfBounds
1104 [ + - ]: 247 : pdfn.writeTxt( pdf[2],
1105 [ - + ][ - - ]: 247 : tk::ctr::PDFInfo{ {"PDF"}, {}, {"ntets"}, 0, 0.0 } );
1106 : :
1107 : 247 : pdfstat_complete(meshid);
1108 [ + - ][ + - ]: 2470 : }
[ + - ][ + + ]
[ - + ][ - + ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ][ - + ]
[ - + ][ + - ]
[ + - ][ + - ]
[ + - ][ + + ]
[ - + ][ - + ]
[ + - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
1109 : :
1110 : : void
1111 [ + - ]: 247 : Transporter::stat()
1112 : : // *****************************************************************************
1113 : : // Echo diagnostics on mesh statistics
1114 : : // *****************************************************************************
1115 : : {
1116 : : tk::Print print;
1117 : :
1118 [ + - ]: 247 : if (++m_nstat == m_nelem.size()) { // stats from all meshes have arrived
1119 : 247 : m_nstat = 0;
1120 : :
1121 : : auto& t = tk::ref_find( m_timer, TimerTag::MESH_PART );
1122 : : print << '\n';
1123 [ + - ]: 741 : print << "Mesh partitioning time: " + std::to_string(t.second) + " sec\n";
1124 : : t = tk::ref_find( m_timer, TimerTag::MESH_DIST );
1125 [ + - ]: 247 : print << "Mesh distribution time: " + std::to_string(t.second) + " sec\n";
1126 : :
1127 : :
1128 [ + + ]: 494 : for (std::size_t i=0; i<m_nelem.size(); ++i) {
1129 [ - + ]: 247 : if (m_nelem.size() > 1) {
1130 [ - - ]: 0 : print.section("Mesh " + std::to_string(i) + " distribution statistics");
1131 : : } else {
1132 [ + - ]: 494 : print.section( "Mesh distribution statistics" );
1133 : : }
1134 : : print <<
1135 : 247 : "min/max/avg edgelength = " +
1136 [ + - ][ + - ]: 741 : std::to_string( m_minstat[i][0] ) + " / " +
[ - + ][ - - ]
1137 [ + - ][ + - ]: 741 : std::to_string( m_maxstat[i][0] ) + " / " +
[ - + ][ - - ]
1138 [ + - ][ - + ]: 494 : std::to_string( m_avgstat[i][0] ) + "\n" +
[ - - ]
1139 [ + - ]: 247 : "min/max/avg V^(1/3) = " +
1140 [ + - ][ + - ]: 741 : std::to_string( m_minstat[i][1] ) + " / " +
[ - + ][ - - ]
1141 [ + - ][ + - ]: 741 : std::to_string( m_maxstat[i][1] ) + " / " +
[ - + ][ - - ]
1142 [ + - ][ - + ]: 494 : std::to_string( m_avgstat[i][1] ) + "\n" +
[ - - ]
1143 [ + - ]: 247 : "min/max/avg nelem = " +
1144 [ + - ][ + - ]: 741 : std::to_string( static_cast<std::size_t>(m_minstat[i][2]) ) + " / " +
[ - + ][ - - ]
1145 [ + - ][ + - ]: 741 : std::to_string( static_cast<std::size_t>(m_maxstat[i][2]) ) + " / " +
[ - + ][ - - ]
1146 [ + - ][ - + ]: 494 : std::to_string( static_cast<std::size_t>(m_avgstat[i][2]) ) + "\n" +
[ - - ]
1147 [ + - ]: 247 : "min/max/avg npoin = " +
1148 [ + - ][ + - ]: 741 : std::to_string( static_cast<std::size_t>(m_minstat[i][3]) ) + " / " +
[ - + ][ - - ]
1149 [ + - ][ + - ]: 741 : std::to_string( static_cast<std::size_t>(m_maxstat[i][3]) ) + " / " +
[ - + ][ - - ]
1150 [ + - ][ - + ]: 494 : std::to_string( static_cast<std::size_t>(m_avgstat[i][3]) ) + "\n" +
[ - - ]
1151 [ + - ]: 247 : "min/max/avg nedge = " +
1152 [ + - ][ + - ]: 741 : std::to_string( static_cast<std::size_t>(m_minstat[i][4]) ) + " / " +
[ - + ][ - - ]
1153 [ + - ][ + - ]: 741 : std::to_string( static_cast<std::size_t>(m_maxstat[i][4]) ) + " / " +
[ - + ][ - - ]
1154 [ + - ][ + - ]: 741 : std::to_string( static_cast<std::size_t>(m_avgstat[i][4]) ) + '\n' +
[ - + ][ - - ]
1155 [ + - ]: 247 : "min/max/avg ncompoin/npoin = " +
1156 [ + - ][ + - ]: 741 : std::to_string( m_minstat[i][5] ) + " / " +
[ - + ][ - - ]
1157 [ + - ][ + - ]: 741 : std::to_string( m_maxstat[i][5] ) + " / " +
[ - + ][ - - ]
1158 [ + - ][ + - ]: 741 : std::to_string( m_avgstat[i][5] ) + '\n';
1159 : : }
1160 : :
1161 : : // Print out time integration header to screen
1162 : 247 : inthead( print );
1163 : :
1164 : 247 : m_progWork.start( print, "Preparing workers", {{ m_nchare[0] }} );
1165 : :
1166 : : // Create "derived-class" workers
1167 [ + + ]: 494 : for (std::size_t i=0; i<m_nelem.size(); ++i) m_sorter[i].createWorkers();
1168 : : }
1169 : 247 : }
1170 : :
1171 : : void
1172 : 247 : Transporter::boxvol( tk::real v, tk::real summeshid )
1173 : : // *****************************************************************************
1174 : : // Reduction target computing total volume of IC box(es)
1175 : : //! \param[in] v Total volume within user-specified IC box(es)
1176 : : //! \param[in] summeshid Mesh id as a real (summed accross the distributed mesh)
1177 : : // *****************************************************************************
1178 : : {
1179 [ + + ]: 247 : auto meshid = tk::cref_find( m_meshid, static_cast<std::size_t>(summeshid) );
1180 [ + + ][ + - ]: 253 : if (v > 0.0) tk::Print() << "IC-box-volume sum: " + std::to_string(v) << '\n';
1181 : :
1182 : : const auto& solver = g_cfg.get< tag::solver >();
1183 [ + + ]: 247 : if (solver == "riecg") {
1184 : 73 : m_riecg[ meshid ].setup( v );
1185 : : }
1186 [ + + ]: 174 : else if (solver == "laxcg") {
1187 : 23 : m_laxcg[ meshid ].setup( v );
1188 : : }
1189 [ + + ]: 151 : else if (solver == "zalcg") {
1190 : 36 : m_zalcg[ meshid ].setup( v );
1191 : : }
1192 [ + + ]: 115 : else if (solver == "kozcg") {
1193 : 44 : m_kozcg[ meshid ].setup( v );
1194 : : }
1195 [ + + ]: 71 : else if (solver == "chocg") {
1196 : 41 : m_chocg[ meshid ].setup( v );
1197 : : }
1198 [ + - ]: 30 : else if (solver == "lohcg") {
1199 : 30 : m_lohcg[ meshid ].setup( v );
1200 : : }
1201 : : else {
1202 [ - - ][ - - ]: 0 : Throw( "Unknown solver: " + solver );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
1203 : : }
1204 : :
1205 : : // Turn on automatic load balancing
1206 [ + - ]: 247 : if (++m_ncom == m_nelem.size()) { // all worker arrays have finished
1207 : 247 : m_ncom = 0;
1208 : : tk::Print print;
1209 : 247 : m_progWork.end( print );
1210 : 247 : tk::CProxy_LBSwitch::ckNew();
1211 : : }
1212 : 247 : }
1213 : :
1214 : : void
1215 : 259 : Transporter::inthead( const tk::Print& print )
1216 : : // *****************************************************************************
1217 : : // Print out time integration header to screen
1218 : : //! \param[in] print Pretty printer object to use for printing
1219 : : // *****************************************************************************
1220 : : {
1221 : 259 : const auto dea = g_cfg.get< tag::deactivate >();
1222 : : const auto solver = g_cfg.get< tag::solver >();
1223 [ + + ]: 259 : const auto pre = solver == "chocg" ? 1 : 0;
1224 : 259 : const auto theta = g_cfg.get< tag::theta >();
1225 : : const auto eps = std::numeric_limits< tk::real >::epsilon();
1226 [ + + ][ + + ]: 259 : const auto mom = solver == "chocg" and theta > eps ? 1 : 0;
1227 : :
1228 [ + - ][ + - ]: 518 : print.section( "Time integration" );
[ + - ][ - - ]
1229 : : print <<
1230 : : "Legend: it - iteration count\n"
1231 : : " t - physics time\n"
1232 : : " dt - physics time step size\n"
1233 : : " ETE - estimated wall-clock time elapsed (h:m:s)\n"
1234 : : " ETA - estimated wall-clock time for accomplishment (h:m:s)\n"
1235 : : " EGT - estimated grind wall-clock time (1e-6sec/timestep)\n"
1236 : : " EGP - estimated grind performance: wall-clock time "
1237 : : "(1e-6sec/DOF/timestep)\n"
1238 : : " flg - status flags, legend:\n"
1239 : : " f - field (volume and surface) output\n"
1240 : : " i - integral output\n"
1241 : : " d - diagnostics output\n"
1242 : : " t - physics time history output\n"
1243 : : " h - h-refinement\n"
1244 : : " l - load balancing\n"
1245 [ + + ]: 259 : " c - checkpoint\n" << (dea ?
1246 [ + + ]: 259 : " e:x/y - x of y work units deactivated\n" : "") << (pre ?
1247 [ + + ]: 259 : " p:it - pressure linear solve iterations\n" : "") << (mom ?
1248 : : " m:it - momentum/transport linear solve iterations\n" : "") <<
1249 : : "\n it t dt ETE ETA EGT"
1250 : : " EGP flg\n"
1251 : : "-----------------------------------------------------------------------"
1252 : : "-----------------\n";
1253 : 259 : }
1254 : :
1255 : : void
1256 [ + - ]: 3201 : Transporter::rhodiagnostics( CkReductionMsg* msg )
1257 : : // *****************************************************************************
1258 : : // Reduction target collecting diagnostics from density-based solvers
1259 : : //! \param[in] msg Serialized diagnostics vector aggregated across all PEs
1260 : : // *****************************************************************************
1261 : : {
1262 : : using namespace diagnostics;
1263 : :
1264 : : std::size_t meshid;
1265 : : std::size_t ncomp;
1266 : : std::vector< std::vector< tk::real > > d;
1267 : :
1268 : : // Deserialize diagnostics vector
1269 : : PUP::fromMem creator( msg->getData() );
1270 : : // cppcheck-suppress uninitvar
1271 : : creator | meshid;
1272 : : creator | ncomp;
1273 : : creator | d;
1274 : : delete msg;
1275 : :
1276 : : // cppcheck-suppress uninitvar
1277 : : // cppcheck-suppress unreadVariable
1278 [ + - ]: 3201 : auto id = std::to_string(meshid);
1279 : :
1280 : : Assert( ncomp > 0, "Number of scalar components must be positive");
1281 : : Assert( d.size() == NUMDIAG, "Diagnostics vector size mismatch" );
1282 : :
1283 : : // cppcheck-suppress unsignedLessThanZero
1284 [ + + ]: 28809 : for (std::size_t i=0; i<d.size(); ++i) {
1285 : : Assert( d[i].size() == ncomp, "Size mismatch at final stage of "
1286 : : "diagnostics aggregation for mesh " + id );
1287 : : }
1288 : :
1289 : : // Allocate storage for those diagnostics that are always computed
1290 [ + - ][ - - ]: 3201 : std::vector< tk::real > diag( ncomp, 0.0 );
1291 : :
1292 : : // Finish computing the L2 norm of conserved variables
1293 [ + + ]: 19480 : for (std::size_t i=0; i<d[L2SOL].size(); ++i) {
1294 : : // cppcheck-suppress uninitvar
1295 : 16279 : diag[i] = sqrt( d[L2SOL][i] / m_meshvol[meshid] );
1296 : : }
1297 : :
1298 : : // Finish computing the L2 norm of the residuals
1299 [ + - ][ - - ]: 3201 : std::vector< tk::real > l2res( d[L2RES].size(), 0.0 );
1300 [ + + ]: 19480 : for (std::size_t i=0; i<d[L2RES].size(); ++i) {
1301 : : // cppcheck-suppress uninitvar
1302 : 16279 : l2res[i] = std::sqrt( d[L2RES][i] / m_meshvol[meshid] );
1303 [ + - ]: 16279 : diag.push_back( l2res[i] );
1304 : : }
1305 : :
1306 : : // Append total energy
1307 [ + - ]: 3201 : diag.push_back( d[TOTALEN][0] );
1308 : :
1309 : : // Finish computing norms of the numerical - analytical solution
1310 [ + - ][ + + ]: 6402 : if (problems::SOL()) {
1311 [ + + ]: 10200 : for (std::size_t i=0; i<d[L2ERR].size(); ++i) {
1312 : : // cppcheck-suppress uninitvar
1313 [ + - ]: 8544 : diag.push_back( std::sqrt( d[L2ERR][i] / m_meshvol[meshid] ) );
1314 : : }
1315 [ + + ]: 10200 : for (std::size_t i=0; i<d[L1ERR].size(); ++i) {
1316 : : // cppcheck-suppress uninitvar
1317 [ + - ][ - - ]: 8544 : diag.push_back( d[L1ERR][i] / m_meshvol[meshid] );
1318 : : }
1319 : : }
1320 : :
1321 : : // Append diagnostics file at selected times
1322 : : auto filename = g_cfg.get< tag::diag >();
1323 [ - + ][ - - ]: 3201 : if (m_nelem.size() > 1) filename += '.' + id;
[ - - ]
1324 : : tk::DiagWriter dw( filename,
1325 : : g_cfg.get< tag::diag_format >(),
1326 : : g_cfg.get< tag::diag_precision >(),
1327 [ + - ]: 3201 : std::ios_base::app );
1328 [ + - ]: 3201 : dw.write( static_cast<uint64_t>(d[ITER][0]), d[TIME][0], d[DT][0], diag );
1329 : :
1330 : : const auto& solver = g_cfg.get< tag::solver >();
1331 [ + + ]: 3201 : if (solver == "riecg") {
1332 : : // cppcheck-suppress uninitvar
1333 [ + - ]: 1613 : m_riecg[ meshid ].evalres( l2res );
1334 : : }
1335 [ + + ]: 1588 : else if (solver == "laxcg") {
1336 : : // cppcheck-suppress uninitvar
1337 [ + - ]: 260 : m_laxcg[ meshid ].evalres( l2res );
1338 : : }
1339 [ + + ]: 1328 : else if (solver == "zalcg") {
1340 : : // cppcheck-suppress uninitvar
1341 [ + - ]: 480 : m_zalcg[ meshid ].evalres( l2res );
1342 : : }
1343 [ + - ]: 848 : else if (solver == "kozcg") {
1344 : : // cppcheck-suppress uninitvar
1345 [ + - ]: 848 : m_kozcg[ meshid ].evalres( l2res );
1346 : : }
1347 : : else {
1348 [ - - ][ - - ]: 0 : Throw( "Unknown solver: " + solver );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
1349 : : }
1350 : 3201 : }
1351 : :
1352 : : void
1353 [ + - ]: 457 : Transporter::prediagnostics( CkReductionMsg* msg )
1354 : : // *****************************************************************************
1355 : : // Reduction target collecting diagnostics from pressure-based solvers
1356 : : //! \param[in] msg Serialized diagnostics vector aggregated across all PEs
1357 : : // *****************************************************************************
1358 : : {
1359 : : using namespace diagnostics;
1360 : :
1361 : : std::size_t meshid;
1362 : : std::size_t ncomp;
1363 : : std::vector< std::vector< tk::real > > d;
1364 : :
1365 : : // Deserialize diagnostics vector
1366 : : PUP::fromMem creator( msg->getData() );
1367 : : // cppcheck-suppress uninitvar
1368 : : creator | meshid;
1369 : : creator | ncomp;
1370 : : creator | d;
1371 : : delete msg;
1372 : :
1373 : : // cppcheck-suppress uninitvar
1374 : : // cppcheck-suppress unreadVariable
1375 [ + - ]: 457 : auto id = std::to_string(meshid);
1376 : :
1377 : : Assert( ncomp > 0, "Number of scalar components must be positive");
1378 : : Assert( d.size() == NUMDIAG, "Diagnostics vector size mismatch" );
1379 : :
1380 : : // cppcheck-suppress unsignedLessThanZero
1381 [ + + ]: 4113 : for (std::size_t i=0; i<d.size(); ++i) {
1382 : : Assert( d[i].size() == ncomp, "Size mismatch at final stage of "
1383 : : "diagnostics aggregation for mesh " + id );
1384 : : }
1385 : :
1386 : : // Allocate storage for those diagnostics that are always computed
1387 [ + - ][ - - ]: 457 : std::vector< tk::real > diag( ncomp, 0.0 );
1388 : :
1389 : : // Finish computing the L2 norm of conserved variables
1390 [ + + ]: 2264 : for (std::size_t i=0; i<d[L2SOL].size(); ++i) {
1391 : : // cppcheck-suppress uninitvar
1392 : 1807 : diag[i] = sqrt( d[L2SOL][i] / m_meshvol[meshid] );
1393 : : }
1394 : :
1395 : : // Finish computing the L2 norm of the residuals
1396 [ + - ][ - - ]: 457 : std::vector< tk::real > l2res( d[L2RES].size(), 0.0 );
1397 [ + + ]: 2264 : for (std::size_t i=0; i<d[L2RES].size(); ++i) {
1398 : : // cppcheck-suppress uninitvar
1399 : 1807 : l2res[i] = std::sqrt( d[L2RES][i] / m_meshvol[meshid] );
1400 [ + - ]: 1807 : diag.push_back( l2res[i] );
1401 : : }
1402 : :
1403 : : // Finish computing norms of the numerical - analytical solution
1404 [ + - ][ + + ]: 914 : if (problems::PRESSURE_SOL()) {
1405 [ + + ]: 14 : for (std::size_t i=0; i<d[L2ERR].size(); ++i) {
1406 : : // cppcheck-suppress uninitvar
1407 [ + - ]: 7 : diag.push_back( std::sqrt( d[L2ERR][i] / m_meshvol[meshid] ) );
1408 : : }
1409 [ + + ]: 14 : for (std::size_t i=0; i<d[L1ERR].size(); ++i) {
1410 : : // cppcheck-suppress uninitvar
1411 [ + - ][ - - ]: 7 : diag.push_back( d[L1ERR][i] / m_meshvol[meshid] );
1412 : : }
1413 : : }
1414 : :
1415 : : // Append diagnostics file at selected times
1416 : : auto filename = g_cfg.get< tag::diag >();
1417 [ - + ][ - - ]: 457 : if (m_nelem.size() > 1) filename += '.' + id;
[ - - ]
1418 : : tk::DiagWriter dw( filename,
1419 : : g_cfg.get< tag::diag_format >(),
1420 : : g_cfg.get< tag::diag_precision >(),
1421 [ + - ]: 457 : std::ios_base::app );
1422 [ + - ]: 457 : dw.write( static_cast<uint64_t>(d[ITER][0]), d[TIME][0], d[DT][0], diag );
1423 : :
1424 : : const auto& solver = g_cfg.get< tag::solver >();
1425 [ + - ]: 457 : if (solver == "chocg") {
1426 : : // cppcheck-suppress uninitvar
1427 [ + - ]: 457 : m_chocg[ meshid ].evalres( l2res );
1428 : : }
1429 : : else {
1430 [ - - ][ - - ]: 0 : Throw( "Unknown solver: " + solver );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
1431 : : }
1432 : 457 : }
1433 : :
1434 : : void
1435 [ + - ]: 364 : Transporter::acdiagnostics( CkReductionMsg* msg )
1436 : : // *****************************************************************************
1437 : : // Reduction target collecting diagnostics from artificial compressibility
1438 : : // solvers
1439 : : //! \param[in] msg Serialized diagnostics vector aggregated across all PEs
1440 : : // *****************************************************************************
1441 : : {
1442 : : using namespace diagnostics;
1443 : :
1444 : : std::size_t meshid;
1445 : : std::size_t ncomp;
1446 : : std::vector< std::vector< tk::real > > d;
1447 : :
1448 : : // Deserialize diagnostics vector
1449 : : PUP::fromMem creator( msg->getData() );
1450 : : // cppcheck-suppress uninitvar
1451 : : creator | meshid;
1452 : : creator | ncomp;
1453 : : creator | d;
1454 : : delete msg;
1455 : :
1456 : : // cppcheck-suppress uninitvar
1457 : : // cppcheck-suppress unreadVariable
1458 [ + - ]: 364 : auto id = std::to_string(meshid);
1459 : :
1460 : : Assert( ncomp > 0, "Number of scalar components must be positive");
1461 : : Assert( d.size() == NUMDIAG, "Diagnostics vector size mismatch" );
1462 : :
1463 : : // cppcheck-suppress unsignedLessThanZero
1464 [ + + ]: 3276 : for (std::size_t i=0; i<d.size(); ++i) {
1465 : : Assert( d[i].size() == ncomp, "Size mismatch at final stage of "
1466 : : "diagnostics aggregation for mesh " + id );
1467 : : }
1468 : :
1469 : : // Allocate storage for those diagnostics that are always computed
1470 [ + - ][ - - ]: 364 : std::vector< tk::real > diag( ncomp, 0.0 );
1471 : :
1472 : : // Finish computing the L2 norm of conserved variables
1473 [ + + ]: 1820 : for (std::size_t i=0; i<d[L2SOL].size(); ++i) {
1474 : : // cppcheck-suppress uninitvar
1475 : 1456 : diag[i] = sqrt( d[L2SOL][i] / m_meshvol[meshid] );
1476 : : }
1477 : :
1478 : : // Finish computing the L2 norm of the residuals
1479 [ + - ][ - - ]: 364 : std::vector< tk::real > l2res( d[L2RES].size(), 0.0 );
1480 [ + + ]: 1820 : for (std::size_t i=0; i<d[L2RES].size(); ++i) {
1481 : : // cppcheck-suppress uninitvar
1482 : 1456 : l2res[i] = std::sqrt( d[L2RES][i] / m_meshvol[meshid] );
1483 [ + - ]: 1456 : diag.push_back( l2res[i] );
1484 : : }
1485 : :
1486 : : // Append diagnostics file at selected times
1487 : : auto filename = g_cfg.get< tag::diag >();
1488 [ - + ][ - - ]: 364 : if (m_nelem.size() > 1) filename += '.' + id;
[ - - ]
1489 : : tk::DiagWriter dw( filename,
1490 : : g_cfg.get< tag::diag_format >(),
1491 : : g_cfg.get< tag::diag_precision >(),
1492 [ + - ]: 364 : std::ios_base::app );
1493 [ + - ]: 364 : dw.write( static_cast<uint64_t>(d[ITER][0]), d[TIME][0], d[DT][0], diag );
1494 : :
1495 : : const auto& solver = g_cfg.get< tag::solver >();
1496 [ + - ]: 364 : if (solver == "lohcg") {
1497 : : // cppcheck-suppress uninitvar
1498 [ + - ]: 364 : m_lohcg[ meshid ].evalres( l2res );
1499 : : }
1500 : : else {
1501 [ - - ][ - - ]: 0 : Throw( "Unknown solver: " + solver );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
1502 : : }
1503 : 364 : }
1504 : :
1505 : : void
1506 [ + - ]: 98 : Transporter::integrals( CkReductionMsg* msg )
1507 : : // *****************************************************************************
1508 : : // Reduction target optionally collecting integrals
1509 : : //! \param[in] msg Serialized integrals aggregated across all PEs
1510 : : // *****************************************************************************
1511 : : {
1512 : : using namespace integrals;
1513 : :
1514 : : // cppcheck-suppress unassignedVariable
1515 : : std::size_t meshid;
1516 : : std::vector< std::map< int, tk::real > > d;
1517 : :
1518 : : // Deserialize integrals vector
1519 : : PUP::fromMem creator( msg->getData() );
1520 : : // cppcheck-suppress uninitvar
1521 : : creator | meshid;
1522 : : creator | d;
1523 : : delete msg;
1524 : :
1525 : : const auto& sidesets_integral = g_cfg.get< tag::integout >();
1526 : : // cppcheck-suppress
1527 [ + - ]: 98 : if (not sidesets_integral.empty()) {
1528 : :
1529 : : Assert( d.size() == NUMINT, "Integrals vector size mismatch" );
1530 : :
1531 : : // Collect integrals for output
1532 : : std::vector< tk::real > ints;
1533 [ + - ][ + + ]: 246 : for (const auto& [s,m] : d[MASS_FLOW_RATE]) ints.push_back( m );
1534 [ + - ][ + + ]: 122 : for (const auto& [s,m] : d[FORCE_X]) ints.push_back( m );
1535 [ + - ][ + + ]: 122 : for (const auto& [s,m] : d[FORCE_Y]) ints.push_back( m );
1536 [ + - ][ + + ]: 122 : for (const auto& [s,m] : d[FORCE_Z]) ints.push_back( m );
1537 : :
1538 : : // Append integrals file at selected times
1539 [ + - ]: 98 : auto filename = g_cfg.get< tag::output >() + ".int";
1540 : : tk::DiagWriter dw( filename,
1541 : : g_cfg.get< tag::integout_format >(),
1542 : : g_cfg.get< tag::integout_precision >(),
1543 [ + - ]: 98 : std::ios_base::app );
1544 : : // cppcheck-suppress containerOutOfBounds
1545 [ + - ]: 294 : dw.write( static_cast<uint64_t>(tk::cref_find( d[ITER], 0 )),
1546 : : // cppcheck-suppress containerOutOfBounds
1547 : : tk::cref_find( d[TIME], 0 ),
1548 : : // cppcheck-suppress containerOutOfBounds
1549 : : tk::cref_find( d[DT], 0 ),
1550 : : ints );
1551 : : }
1552 : :
1553 : : const auto& solver = g_cfg.get< tag::solver >();
1554 [ + + ]: 98 : if (solver == "riecg") {
1555 : : // cppcheck-suppress uninitvar
1556 [ + - ]: 74 : m_riecg[ meshid ].step();
1557 : : }
1558 [ - + ]: 24 : else if (solver == "laxcg") {
1559 : : // cppcheck-suppress uninitvar
1560 [ - - ]: 0 : m_laxcg[ meshid ].step();
1561 : : }
1562 [ - + ]: 24 : else if (solver == "zalcg") {
1563 : : // cppcheck-suppress uninitvar
1564 [ - - ]: 0 : m_zalcg[ meshid ].step();
1565 : : }
1566 [ - + ]: 24 : else if (solver == "kozcg") {
1567 : : // cppcheck-suppress uninitvar
1568 [ - - ]: 0 : m_kozcg[ meshid ].step();
1569 : : }
1570 [ + + ]: 24 : else if (solver == "chocg") {
1571 : : // cppcheck-suppress uninitvar
1572 [ + - ]: 20 : m_chocg[ meshid ].step();
1573 : : }
1574 [ + - ]: 4 : else if (solver == "lohcg") {
1575 : : // cppcheck-suppress uninitvar
1576 [ + - ]: 4 : m_lohcg[ meshid ].step();
1577 : : }
1578 : : else
1579 [ - - ][ - - ]: 0 : Throw( "Unknown solver: " + solver );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
1580 : 98 : }
1581 : :
1582 : : void
1583 [ + + ]: 271 : Transporter::resume()
1584 : : // *****************************************************************************
1585 : : // Resume execution from checkpoint/restart files
1586 : : //! \details This is invoked by Charm++ after the checkpoint is done, as well as
1587 : : //! when the restart (returning from a checkpoint) is complete
1588 : : // *****************************************************************************
1589 : : {
1590 [ + + ]: 271 : if (std::any_of(begin(m_finished), end(m_finished), [](auto f){return !f;})) {
1591 : :
1592 : : // If just restarted from a checkpoint, Main( CkMigrateMessage* msg ) has
1593 : : // increased g_nrestart, but only on PE 0, so broadcast.
1594 : :
1595 : : const auto& solver = g_cfg.get< tag::solver >();
1596 [ + + ]: 12 : if (solver == "riecg") {
1597 [ + + ]: 4 : for (std::size_t i=0; i<m_nelem.size(); ++i) {
1598 : 2 : m_riecg[i].evalLB( g_nrestart );
1599 : : }
1600 : : }
1601 [ + + ]: 10 : else if (solver == "laxcg") {
1602 [ + + ]: 4 : for (std::size_t i=0; i<m_nelem.size(); ++i) {
1603 : 2 : m_laxcg[i].evalLB( g_nrestart );
1604 : : }
1605 : : }
1606 [ + + ]: 8 : else if (solver == "zalcg") {
1607 [ + + ]: 4 : for (std::size_t i=0; i<m_nelem.size(); ++i) {
1608 : 2 : m_zalcg[i].evalLB( g_nrestart );
1609 : : }
1610 : : }
1611 [ + + ]: 6 : else if ( solver == "kozcg") {
1612 [ + + ]: 4 : for (std::size_t i=0; i<m_nelem.size(); ++i) {
1613 : 2 : m_kozcg[i].evalLB( g_nrestart );
1614 : : }
1615 : : }
1616 [ + + ]: 4 : else if ( solver == "chocg") {
1617 [ + + ]: 4 : for (std::size_t i=0; i<m_nelem.size(); ++i) {
1618 : 2 : m_chocg[i].evalLB( g_nrestart );
1619 : : }
1620 : : }
1621 [ + - ]: 2 : else if ( solver == "lohcg") {
1622 [ + + ]: 4 : for (std::size_t i=0; i<m_nelem.size(); ++i) {
1623 : 2 : m_lohcg[i].evalLB( g_nrestart );
1624 : : }
1625 : : }
1626 : : else {
1627 [ - - ][ - - ]: 0 : Throw( "Unknown solver: " + solver );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
1628 : : }
1629 : :
1630 : :
1631 : : } else {
1632 : :
1633 : 259 : mainProxy.finalize();
1634 : :
1635 : : }
1636 : 271 : }
1637 : :
1638 : : void
1639 [ + - ]: 259 : Transporter::checkpoint( std::size_t finished, std::size_t meshid )
1640 : : // *****************************************************************************
1641 : : // Save checkpoint/restart files
1642 : : //! \param[in] finished Nonzero if finished with time stepping
1643 : : //! \param[in] meshid Mesh id
1644 : : // *****************************************************************************
1645 : : {
1646 : 259 : m_finished[meshid] = finished;
1647 : :
1648 [ + - ]: 259 : if (++m_nchk == m_nelem.size()) { // all worker arrays have checkpointed
1649 : 259 : m_nchk = 0;
1650 [ + + ]: 259 : if (not g_cfg.get< tag::benchmark >()) {
1651 : : const auto& ckptdir = g_cfg.get< tag::checkpoint >();
1652 : : CkCallback res( CkIndex_Transporter::resume(), thisProxy );
1653 [ + - ]: 151 : CkStartCheckpoint( ckptdir.c_str(), res );
1654 : : //CkStartMemCheckpoint( res );
1655 : : } else {
1656 : 108 : resume();
1657 : : }
1658 : : }
1659 : 259 : }
1660 : :
1661 : : void
1662 : 259 : Transporter::finish( std::size_t meshid )
1663 : : // *****************************************************************************
1664 : : // Normal finish of time stepping
1665 : : //! \param[in] meshid Mesh id
1666 : : // *****************************************************************************
1667 : : {
1668 : 259 : checkpoint( /* finished = */ 1, meshid );
1669 : 259 : }
1670 : :
1671 : : #include "NoWarning/transporter.def.h"
|
