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