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