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