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