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 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Transporter drives the time integration of transport equations
9 : : \details Transporter drives the time integration of transport equations.
10 : : The implementation uses the Charm++ runtime system and is fully asynchronous,
11 : : overlapping computation, communication as well as I/O. The algorithm
12 : : utilizes the structured dagger (SDAG) Charm++ functionality. The high-level
13 : : overview of the algorithm structure and how it interfaces with Charm++ is
14 : : discussed in the Charm++ interface file src/Inciter/transporter.ci.
15 : : */
16 : : // *****************************************************************************
17 : :
18 : : #include <string>
19 : : #include <iostream>
20 : : #include <cstddef>
21 : : #include <unordered_set>
22 : : #include <limits>
23 : : #include <cmath>
24 : :
25 : : #include <brigand/algorithms/for_each.hpp>
26 : :
27 : : #include "Macro.hpp"
28 : : #include "Transporter.hpp"
29 : : #include "Fields.hpp"
30 : : #include "PDEStack.hpp"
31 : : #include "UniPDF.hpp"
32 : : #include "PDFWriter.hpp"
33 : : #include "ContainerUtil.hpp"
34 : : #include "LoadDistributor.hpp"
35 : : #include "MeshReader.hpp"
36 : : #include "Inciter/Types.hpp"
37 : : #include "Inciter/InputDeck/InputDeck.hpp"
38 : : #include "NodeDiagnostics.hpp"
39 : : #include "ElemDiagnostics.hpp"
40 : : #include "DiagWriter.hpp"
41 : : #include "Callback.hpp"
42 : : #include "CartesianProduct.hpp"
43 : :
44 : : #include "NoWarning/inciter.decl.h"
45 : : #include "NoWarning/partitioner.decl.h"
46 : :
47 : : extern CProxy_Main mainProxy;
48 : :
49 : : namespace inciter {
50 : :
51 : : extern ctr::InputDeck g_inputdeck;
52 : : extern ctr::InputDeck g_inputdeck_defaults;
53 : : extern std::vector< CGPDE > g_cgpde;
54 : : extern std::vector< DGPDE > g_dgpde;
55 : :
56 : : }
57 : :
58 : : using inciter::Transporter;
59 : :
60 : 233 : Transporter::Transporter() :
61 : : m_input( input() ),
62 : : m_nchare( m_input.size() ),
63 : : m_meshid(),
64 : 0 : m_ncit( m_nchare.size(), 0 ),
65 : : m_nload( 0 ),
66 : : m_npart( 0 ),
67 : : m_nstat( 0 ),
68 : : m_ndisc( 0 ),
69 : : m_nchk( 0 ),
70 : : m_ncom( 0 ),
71 : 0 : m_nt0refit( m_nchare.size(), 0 ),
72 : 0 : m_ndtrefit( m_nchare.size(), 0 ),
73 : 0 : m_noutrefit( m_nchare.size(), 0 ),
74 : 0 : m_noutderefit( m_nchare.size(), 0 ),
75 : : m_scheme(),
76 : : m_partitioner(),
77 : : m_refiner(),
78 : : m_meshwriter(),
79 : : m_sorter(),
80 : : m_nelem( m_nchare.size() ),
81 : : m_npoin(),
82 : 0 : m_finished( m_nchare.size(), 0 ),
83 : : m_meshvol( m_nchare.size() ),
84 : : m_minstat( m_nchare.size() ),
85 : : m_maxstat( m_nchare.size() ),
86 : : m_avgstat( m_nchare.size() ),
87 : : m_timer(),
88 : 466 : m_progMesh( g_inputdeck.get< tag::cmd, tag::feedback >(),
89 : : ProgMeshPrefix, ProgMeshLegend ),
90 : 233 : m_progWork( g_inputdeck.get< tag::cmd, tag::feedback >(),
91 [ + - ][ + - ]: 1398 : ProgWorkPrefix, ProgWorkLegend )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ][ + - ]
92 : : // *****************************************************************************
93 : : // Constructor
94 : : // *****************************************************************************
95 : : {
96 : : // Echo configuration to screen
97 [ + - ][ + - ]: 233 : info( printer() );
98 : :
99 : 233 : const auto nstep = g_inputdeck.get< tag::discr, tag::nstep >();
100 : 233 : const auto t0 = g_inputdeck.get< tag::discr, tag::t0 >();
101 : 233 : const auto term = g_inputdeck.get< tag::discr, tag::term >();
102 : 233 : const auto constdt = g_inputdeck.get< tag::discr, tag::dt >();
103 : :
104 : : // If the desired max number of time steps is larger than zero, and the
105 : : // termination time is larger than the initial time, and the constant time
106 : : // step size (if that is used) is smaller than the duration of the time to be
107 : : // simulated, we have work to do, otherwise, finish right away. If a constant
108 : : // dt is not used, that part of the logic is always true as the default
109 : : // constdt is zero, see inciter::ctr::InputDeck::InputDeck().
110 [ + - ][ + - ]: 233 : if ( nstep != 0 && term > t0 && constdt < term-t0 ) {
[ + - ]
111 : :
112 : : // Enable SDAG waits for collecting mesh statistics
113 [ + - ]: 233 : thisProxy.wait4stat();
114 : :
115 : : // Configure and write diagnostics file header
116 [ + - ]: 233 : diagHeader();
117 : :
118 : : // Create mesh partitioner AND boundary condition object group
119 [ + - ]: 233 : createPartitioner();
120 : :
121 [ - - ]: 0 : } else finish(); // stop if no time stepping requested
122 : 233 : }
123 : :
124 : 0 : Transporter::Transporter( CkMigrateMessage* m ) :
125 : : CBase_Transporter( m ),
126 : 0 : m_progMesh( g_inputdeck.get< tag::cmd, tag::feedback >(),
127 : : ProgMeshPrefix, ProgMeshLegend ),
128 : 0 : m_progWork( g_inputdeck.get< tag::cmd, tag::feedback >(),
129 [ - - ][ - - ]: 0 : ProgWorkPrefix, ProgWorkLegend )
[ - - ]
130 : : // *****************************************************************************
131 : : // Migrate constructor: returning from a checkpoint
132 : : //! \param[in] m Charm++ migrate message
133 : : // *****************************************************************************
134 : : {
135 [ - - ]: 0 : auto print = printer();
136 [ - - ][ - - ]: 0 : print.diag( "Restarted from checkpoint" );
137 [ - - ]: 0 : info( print );
138 [ - - ]: 0 : inthead( print );
139 : 0 : }
140 : :
141 : : std::vector< std::string >
142 : 233 : Transporter::input()
143 : : // *****************************************************************************
144 : : // Generate list of input mesh filenames configured by the user
145 : : //! \return List of input mesh filenames configured by the user
146 : : //! \details If the input file is given on the command line, a single solver
147 : : //! will be instantiated on the single mesh, solving potentially multiple
148 : : //! systems of (potentially coupled) equations. If the input file is not given
149 : : //! on the command line, the mesh files are expected to be configured in the
150 : : //! control/input file, associating a potentially different mesh to each
151 : : //! solver. Both configurations allow the solution of coupled systems, but the
152 : : //! first one solves all equations on the same mesh, while the latter can
153 : : //! couple solutions computed on multiple different meshes.
154 : : // *****************************************************************************
155 : : {
156 : : // Query input mesh filename specified on the command line
157 : 233 : const auto& cmdinput = g_inputdeck.get< tag::cmd, tag::io, tag::input >();
158 : :
159 : : // Extract mesh filenames specified in the control file (assigned to solvers)
160 [ + - ]: 466 : auto ctrinput = g_inputdeck.mesh();
161 : :
162 [ + + ][ - + ]: 233 : ErrChk( !cmdinput.empty() || !ctrinput.empty(),
[ - - ][ - - ]
[ - - ]
163 : : "Either a single input mesh must be given on the command line or multiple "
164 : : "meshes must be configured in the control file." );
165 : :
166 : : // Prepend control file path to mesh filenames in given in control file
167 [ + + ]: 233 : if (!ctrinput.empty()) {
168 : 10 : const auto& ctr = g_inputdeck.get< tag::cmd, tag::io, tag::control >();
169 [ + - ]: 20 : auto path = ctr.substr( 0, ctr.find_last_of("/")+1 );
170 [ + + ][ + - ]: 20 : for (auto& f : ctrinput) f = path + f;
171 : : }
172 : :
173 [ + + ][ + - ]: 458 : if (cmdinput.empty()) return ctrinput; else return { cmdinput };
[ + - ][ + + ]
[ - - ]
174 : : }
175 : :
176 : : void
177 : 233 : Transporter::info( const InciterPrint& print )
178 : : // *****************************************************************************
179 : : // Echo configuration to screen
180 : : //! \param[in] print Pretty printer object to use for printing
181 : : // *****************************************************************************
182 : : {
183 [ + - ][ + - ]: 233 : print.part( "Factory" );
184 : :
185 : : // Print out info data layout
186 [ + - ][ + - ]: 233 : print.list( "Unknowns data layout (CMake: FIELD_DATA_LAYOUT)",
187 [ + - ][ + - ]: 699 : std::list< std::string >{ tk::Fields::layout() } );
[ + + ][ - - ]
188 : :
189 : : // Re-create partial differential equations stack for output
190 [ + - ]: 466 : PDEStack stack;
191 : :
192 : : // Print out information on PDE factories
193 [ + - ][ + - ]: 466 : print.eqlist( "Registered PDEs using continuous Galerkin (CG) methods",
194 : 233 : stack.cgfactory(), stack.cgntypes() );
195 [ + - ][ + - ]: 466 : print.eqlist( "Registered PDEs using discontinuous Galerkin (DG) methods",
196 : 233 : stack.dgfactory(), stack.dgntypes() );
197 [ + - ]: 233 : print.endpart();
198 : :
199 : : // Print out information on problem
200 [ + - ][ + - ]: 233 : print.part( "Problem" );
201 : :
202 : : // Print out info on problem title
203 [ + - ]: 233 : if ( !g_inputdeck.get< tag::title >().empty() )
204 [ + - ]: 233 : print.title( g_inputdeck.get< tag::title >() );
205 : :
206 : 233 : const auto nstep = g_inputdeck.get< tag::discr, tag::nstep >();
207 : 233 : const auto t0 = g_inputdeck.get< tag::discr, tag::t0 >();
208 : 233 : const auto term = g_inputdeck.get< tag::discr, tag::term >();
209 : 233 : const auto constdt = g_inputdeck.get< tag::discr, tag::dt >();
210 : 233 : const auto cfl = g_inputdeck.get< tag::discr, tag::cfl >();
211 : 233 : const auto scheme = g_inputdeck.get< tag::discr, tag::scheme >();
212 : :
213 : : // Print discretization parameters
214 [ + - ][ + - ]: 233 : print.section( "Discretization parameters" );
215 [ + - ]: 233 : print.Item< ctr::Scheme, tag::discr, tag::scheme >();
216 : :
217 [ + + ]: 233 : if (scheme == ctr::SchemeType::DiagCG) {
218 : 83 : auto fct = g_inputdeck.get< tag::discr, tag::fct >();
219 [ + - ][ + - ]: 83 : print.item( "Flux-corrected transport (FCT)", fct );
220 [ + + ]: 83 : if (fct) {
221 [ + - ][ + - ]: 82 : print.item( "FCT mass diffusion coeff",
222 : 82 : g_inputdeck.get< tag::discr, tag::ctau >() );
223 [ + - ][ + - ]: 82 : print.item( "FCT small number",
224 : 82 : g_inputdeck.get< tag::discr, tag::fcteps >() );
225 [ + - ][ + - ]: 82 : print.item( "Clipping FCT",
226 : 82 : g_inputdeck.get< tag::discr, tag::fctclip >() );
227 : : }
228 [ + + ][ + + ]: 150 : } else if (scheme == ctr::SchemeType::DG ||
229 [ + + ][ + + ]: 92 : scheme == ctr::SchemeType::P0P1 || scheme == ctr::SchemeType::DGP1 ||
230 [ + + ]: 66 : scheme == ctr::SchemeType::DGP2 || scheme == ctr::SchemeType::PDG)
231 : : {
232 [ + - ]: 96 : print.Item< ctr::Limiter, tag::discr, tag::limiter >();
233 : : }
234 [ + - ][ + - ]: 233 : print.item( "PE-locality mesh reordering",
235 : 233 : g_inputdeck.get< tag::discr, tag::pelocal_reorder >() );
236 [ + - ][ + - ]: 233 : print.item( "Operator-access mesh reordering",
237 : 233 : g_inputdeck.get< tag::discr, tag::operator_reorder >() );
238 : 233 : auto steady = g_inputdeck.get< tag::discr, tag::steady_state >();
239 [ + - ][ + - ]: 233 : print.item( "Local time stepping", steady );
240 [ + + ]: 233 : if (steady) {
241 [ + - ][ + - ]: 2 : print.item( "L2-norm residual convergence criterion",
242 : 2 : g_inputdeck.get< tag::discr, tag::residual >() );
243 [ + - ][ + - ]: 2 : print.item( "Convergence criterion component index",
244 : 2 : g_inputdeck.get< tag::discr, tag::rescomp >() );
245 : : }
246 [ + - ][ + - ]: 233 : print.item( "Number of time steps", nstep );
247 [ + - ][ + - ]: 233 : print.item( "Start time", t0 );
248 [ + - ][ + - ]: 233 : print.item( "Terminate time", term );
249 : :
250 [ + + ]: 466 : if (std::abs(constdt - g_inputdeck_defaults.get< tag::discr, tag::dt >()) >
251 : 233 : std::numeric_limits< tk::real >::epsilon())
252 [ + - ][ + - ]: 117 : print.item( "Constant time step size", constdt );
253 [ + - ]: 232 : else if (std::abs(cfl - g_inputdeck_defaults.get< tag::discr, tag::cfl >()) >
254 : 116 : std::numeric_limits< tk::real >::epsilon())
255 : : {
256 [ + - ][ + - ]: 116 : print.item( "CFL coefficient", cfl );
257 : : }
258 : :
259 : : // Print out info on settings of selected partial differential equations
260 [ + - ][ + - ]: 233 : print.pdes( "Partial differential equations integrated", stack.info() );
[ + - ]
261 : :
262 : : // Print out solver coupling info
263 [ + - ]: 233 : print.couple( g_inputdeck.get< tag::couple, tag::transfer >(),
264 [ + - ]: 466 : g_inputdeck.depvar() );
265 : :
266 : : // Print out adaptive polynomial refinement configuration
267 [ + + ]: 233 : if (scheme == ctr::SchemeType::PDG) {
268 [ + - ][ + - ]: 12 : print.section( "Polynomial refinement (p-ref)" );
269 [ + - ][ + - ]: 12 : print.item( "p-refinement",
270 : 12 : g_inputdeck.get< tag::pref, tag::pref >() );
271 [ + - ]: 12 : print.Item< ctr::PrefIndicator, tag::pref, tag::indicator >();
272 [ + - ][ + - ]: 12 : print.item( "Max degrees of freedom",
273 : 12 : g_inputdeck.get< tag::pref, tag::ndofmax >() );
274 [ + - ][ + - ]: 12 : print.item( "Tolerance",
275 : 12 : g_inputdeck.get< tag::pref, tag::tolref >() );
276 : : }
277 : :
278 : : // Print out adaptive mesh refinement configuration
279 : 233 : const auto amr = g_inputdeck.get< tag::amr, tag::amr >();
280 [ + + ]: 233 : if (amr) {
281 [ + - ][ + - ]: 42 : print.section( "Mesh refinement (h-ref)" );
282 [ + - ]: 42 : print.refvar( g_inputdeck.get< tag::amr, tag::refvar >(),
283 : 42 : g_inputdeck.get< tag::amr, tag::id >() );
284 [ + - ]: 42 : print.Item< ctr::AMRError, tag::amr, tag::error >();
285 : 42 : auto t0ref = g_inputdeck.get< tag::amr, tag::t0ref >();
286 [ + - ][ + - ]: 42 : print.item( "Refinement at t<0 (t0ref)", t0ref );
287 [ + + ]: 42 : if (t0ref) {
288 : 25 : const auto& initref = g_inputdeck.get< tag::amr, tag::init >();
289 [ + - ][ + - ]: 25 : print.item( "Initial refinement steps", initref.size() );
290 [ + - ]: 25 : print.edgeref( g_inputdeck.get< tag::amr, tag::edge >() );
291 : :
292 : 25 : auto eps =
293 : : std::numeric_limits< kw::amr_xminus::info::expect::type >::epsilon();
294 : :
295 : 25 : auto xminus = g_inputdeck.get< tag::amr, tag::xminus >();
296 : 25 : auto xminus_default = g_inputdeck_defaults.get< tag::amr, tag::xminus >();
297 [ + + ]: 25 : if (std::abs( xminus - xminus_default ) > eps)
298 [ + - ][ + - ]: 1 : print.item( "Initial refinement x-", xminus );
299 : 25 : auto xplus = g_inputdeck.get< tag::amr, tag::xplus >();
300 : 25 : auto xplus_default = g_inputdeck_defaults.get< tag::amr, tag::xplus >();
301 [ - + ]: 25 : if (std::abs( xplus - xplus_default ) > eps)
302 [ - - ][ - - ]: 0 : print.item( "Initial refinement x+", xplus );
303 : :
304 : 25 : auto yminus = g_inputdeck.get< tag::amr, tag::yminus >();
305 : 25 : auto yminus_default = g_inputdeck_defaults.get< tag::amr, tag::yminus >();
306 [ - + ]: 25 : if (std::abs( yminus - yminus_default ) > eps)
307 [ - - ][ - - ]: 0 : print.item( "Initial refinement y-", yminus );
308 : 25 : auto yplus = g_inputdeck.get< tag::amr, tag::yplus >();
309 : 25 : auto yplus_default = g_inputdeck_defaults.get< tag::amr, tag::yplus >();
310 [ - + ]: 25 : if (std::abs( yplus - yplus_default ) > eps)
311 [ - - ][ - - ]: 0 : print.item( "Initial refinement y+", yplus );
312 : :
313 : 25 : auto zminus = g_inputdeck.get< tag::amr, tag::zminus >();
314 : 25 : auto zminus_default = g_inputdeck_defaults.get< tag::amr, tag::zminus >();
315 [ - + ]: 25 : if (std::abs( zminus - zminus_default ) > eps)
316 [ - - ][ - - ]: 0 : print.item( "Initial refinement z-", zminus );
317 : 25 : auto zplus = g_inputdeck.get< tag::amr, tag::zplus >();
318 : 25 : auto zplus_default = g_inputdeck_defaults.get< tag::amr, tag::zplus >();
319 [ - + ]: 25 : if (std::abs( zplus - zplus_default ) > eps)
320 [ - - ][ - - ]: 0 : print.item( "Initial refinement z+", zplus );
321 : : }
322 : 42 : auto dtref = g_inputdeck.get< tag::amr, tag::dtref >();
323 [ + - ][ + - ]: 42 : print.item( "Refinement at t>0 (dtref)", dtref );
324 [ + + ]: 42 : if (dtref) {
325 : 17 : auto dtfreq = g_inputdeck.get< tag::amr, tag::dtfreq >();
326 [ + - ][ + - ]: 17 : print.item( "Mesh refinement frequency, t>0", dtfreq );
327 [ + - ][ + - ]: 17 : print.item( "Uniform-only mesh refinement, t>0",
328 : 17 : g_inputdeck.get< tag::amr, tag::dtref_uniform >() );
329 : : }
330 [ + - ][ + - ]: 42 : print.item( "Refinement tolerance",
331 : 42 : g_inputdeck.get< tag::amr, tag::tolref >() );
332 [ + - ][ + - ]: 42 : print.item( "De-refinement tolerance",
333 : 42 : g_inputdeck.get< tag::amr, tag::tolderef >() );
334 : : }
335 : :
336 : : // Print out ALE configuration
337 : 233 : const auto ale = g_inputdeck.get< tag::ale, tag::ale >();
338 [ + + ]: 233 : if (ale) {
339 [ + - ][ + - ]: 10 : print.section( "Arbitrary Lagrangian-Eulerian (ALE) mesh motion" );
340 : 10 : auto dvcfl = g_inputdeck.get< tag::ale, tag::dvcfl >();
341 [ + - ][ + - ]: 10 : print.item( "Volume-change CFL coefficient", dvcfl );
342 [ + - ]: 10 : print.Item< ctr::MeshVelocity, tag::ale, tag::meshvelocity >();
343 [ + - ]: 10 : print.Item< ctr::MeshVelocitySmoother, tag::ale, tag::smoother >();
344 [ + - ][ + - ]: 10 : print.item( "Mesh motion dimensions", tk::parameters(
[ + - ]
345 : 10 : g_inputdeck.get< tag::ale, tag::mesh_motion >() ) );
346 : 10 : const auto& meshforce = g_inputdeck.get< tag::ale, tag::meshforce >();
347 [ + - ][ + - ]: 10 : print.item( "Mesh velocity force coefficients", tk::parameters(meshforce) );
[ + - ]
348 [ + - ][ + - ]: 10 : print.item( "Vorticity multiplier",
349 : 10 : g_inputdeck.get< tag::ale, tag::vortmult >() );
350 [ + - ][ + - ]: 10 : print.item( "Mesh velocity linear solver tolerance",
351 : 10 : g_inputdeck.get< tag::ale, tag::tolerance >() );
352 [ + - ][ + - ]: 10 : print.item( "Mesh velocity linear solver maxit",
353 : 10 : g_inputdeck.get< tag::ale, tag::maxit >() );
354 : 10 : const auto& dir = g_inputdeck.get< tag::ale, tag::bcdir >();
355 [ + + ]: 10 : if (not dir.empty())
356 [ + - ][ + - ]: 5 : print.item( "Mesh velocity Dirichlet BC sideset(s)",
357 [ + - ]: 10 : tk::parameters( dir ) );
358 : 10 : const auto& sym = g_inputdeck.get< tag::ale, tag::bcsym >();
359 [ + + ]: 10 : if (not sym.empty())
360 [ + - ][ + - ]: 3 : print.item( "Mesh velocity symmetry BC sideset(s)",
361 [ + - ]: 6 : tk::parameters( sym ) );
362 : 10 : std::size_t i = 1;
363 [ + + ]: 12 : for (const auto& m : g_inputdeck.get< tag::ale, tag::move >()) {
364 [ + - ]: 2 : tk::ctr::UserTable opt;
365 [ + - ][ + - ]: 2 : print.item( opt.group() + ' ' + std::to_string(i) + " interpreted as",
[ + - ][ + - ]
366 [ + - ][ + - ]: 2 : opt.name( m.get< tag::fntype >() ) );
367 : 2 : const auto& s = m.get< tag::sideset >();
368 [ + - ]: 2 : if (not s.empty())
369 [ + - ][ + - ]: 2 : print.item( "Moving sideset(s) with table " + std::to_string(i),
[ + - ]
370 [ + - ]: 4 : tk::parameters(s));
371 : 2 : ++i;
372 : : }
373 : : }
374 : :
375 : : // Print I/O filenames
376 [ + - ][ + - ]: 233 : print.section( "Input/Output filenames and directories" );
377 [ + - ][ + - ]: 233 : print.item( "Input mesh(es)", tk::parameters( m_input ) );
[ + - ]
378 : 233 : const auto& of = g_inputdeck.get< tag::cmd, tag::io, tag::output >();
379 [ + - ][ + - ]: 233 : print.item( "Volume field output file(s)",
380 [ + - ]: 466 : of + ".e-s.<meshid>.<numchares>.<chareid>" );
381 [ + - ][ + - ]: 233 : print.item( "Surface field output file(s)",
382 [ + - ]: 466 : of + "-surf.<surfid>.e-s.<meshid>.<numchares>.<chareid>" );
383 [ + - ][ + - ]: 233 : print.item( "History output file(s)", of + ".hist.{pointid}" );
[ + - ]
384 [ + - ][ + - ]: 233 : print.item( "Diagnostics file",
385 : 233 : g_inputdeck.get< tag::cmd, tag::io, tag::diag >() );
386 [ + - ][ + - ]: 233 : print.item( "Checkpoint/restart directory",
387 [ + - ]: 466 : g_inputdeck.get< tag::cmd, tag::io, tag::restart >() + '/' );
388 : :
389 : : // Print output intervals
390 [ + - ][ + - ]: 233 : print.section( "Output intervals (in units of iteration count)" );
391 [ + - ][ + - ]: 233 : print.item( "TTY", g_inputdeck.get< tag::output, tag::iter, tag::tty>() );
392 [ + - ][ + - ]: 233 : print.item( "Field and surface",
393 : 233 : g_inputdeck.get< tag::output, tag::iter, tag::field >() );
394 [ + - ][ + - ]: 233 : print.item( "History",
395 : 233 : g_inputdeck.get< tag::output, tag::iter, tag::history >() );
396 [ + - ][ + - ]: 233 : print.item( "Diagnostics",
397 : 233 : g_inputdeck.get< tag::output, tag::iter, tag::diag >() );
398 [ + - ][ + - ]: 233 : print.item( "Checkpoint/restart",
399 : 233 : g_inputdeck.get< tag::cmd, tag::rsfreq >() );
400 : 233 : auto tf = g_inputdeck.get< tag::output, tag::time, tag::field >();
401 : 233 : auto th = g_inputdeck.get< tag::output, tag::time, tag::history >();
402 [ + + ][ + + ]: 233 : if (tf>0.0 || th>0.0) {
403 [ + - ][ + - ]: 4 : print.section( "Output intervals (in units of physics time)" );
404 [ + + ][ + - ]: 4 : if (tf > 0.0) print.item( "Field and surface", tf );
[ + - ]
405 [ + + ][ + - ]: 4 : if (th > 0.0) print.item( "History", th );
[ + - ]
406 : : }
407 : 233 : const auto& rf = g_inputdeck.get< tag::output, tag::range, tag::field >();
408 : 233 : const auto& rh = g_inputdeck.get< tag::output, tag::range, tag::history >();
409 [ + - ][ + + ]: 233 : if (not rf.empty() or not rh.empty()) {
[ + + ]
410 [ + - ][ + - ]: 1 : print.section( "Output time ranges (in units of physics time)" );
411 [ - + ]: 1 : for (const auto& r : rf)
412 [ - - ][ - - ]: 0 : print.item("Field output { mintime, maxtime, dt }", tk::parameters(r));
[ - - ]
413 [ + + ]: 3 : for (const auto& r : rh)
414 [ + - ][ + - ]: 2 : print.item("History output { mintime, maxtime, dt }", tk::parameters(r));
[ + - ]
415 : : }
416 : :
417 : : // Print output variables: fields and surfaces
418 [ + - ]: 466 : const auto nodeoutvars = g_inputdeck.outvars( tk::Centering::NODE );
419 [ + - ]: 466 : const auto elemoutvars = g_inputdeck.outvars( tk::Centering::ELEM );
420 [ + - ]: 466 : const auto aliases = g_inputdeck.outvar_aliases();
421 [ + - ]: 466 : const auto outsets = g_inputdeck.outsets();
422 [ + + ][ + + ]: 233 : if (!nodeoutvars.empty() || !elemoutvars.empty() || !outsets.empty())
[ - + ][ + + ]
423 [ + - ][ + - ]: 165 : print.section( "Output fields" );
424 [ + + ]: 233 : if (!nodeoutvars.empty())
425 [ + - ][ + - ]: 88 : print.item( "Node field(s)", tk::parameters(nodeoutvars) );
[ + - ]
426 [ + + ]: 233 : if (!elemoutvars.empty())
427 [ + - ][ + - ]: 77 : print.item( "Elem field(s)", tk::parameters(elemoutvars) );
[ + - ]
428 [ + + ]: 233 : if (!aliases.empty())
429 [ + - ][ + - ]: 154 : print.item( "Alias(es)", tk::parameters(aliases) );
[ + - ]
430 [ + + ]: 233 : if (!outsets.empty())
431 [ + - ][ + - ]: 4 : print.item( "Surface side set(s)", tk::parameters(outsets) );
[ + - ]
432 : :
433 : : // Print output variables: history
434 : 233 : const auto& pt = g_inputdeck.get< tag::history, tag::point >();
435 : 233 : const auto& id = g_inputdeck.get< tag::history, tag::id >();
436 [ + + ]: 233 : if (!pt.empty()) {
437 [ + - ][ + - ]: 7 : print.section( "Output time history" );
438 [ + + ]: 21 : for (std::size_t p=0; p<pt.size(); ++p) {
439 [ + - ]: 14 : std::stringstream ss;
440 : 14 : auto prec = g_inputdeck.get< tag::prec, tag::history >();
441 [ + - ]: 14 : ss << std::setprecision( static_cast<int>(prec) );
442 [ + - ][ + - ]: 14 : ss << of << ".hist." << id[p];
[ + - ]
443 [ + - ][ + - ]: 14 : print.longitem( "At point " + id[p] + ' ' + tk::parameters(pt[p]),
[ + - ][ + - ]
[ + - ]
444 [ + - ]: 28 : ss.str() );
445 : : }
446 : : }
447 : :
448 [ + - ]: 233 : print.endsubsection();
449 : 233 : }
450 : :
451 : : bool
452 : 299 : Transporter::matchBCs( std::map< int, std::vector< std::size_t > >& bnd )
453 : : // *****************************************************************************
454 : : // Verify that side sets specified in the control file exist in mesh file
455 : : //! \details This function does two things: (1) it verifies that the side
456 : : //! sets used in the input file (either to which boundary conditions (BC)
457 : : //! are assigned or listed as field output by the user in the
458 : : //! input file) all exist among the side sets read from the input mesh
459 : : //! file and errors out if at least one does not, and (2) it matches the
460 : : //! side set ids at which the user has configured BCs (or listed as an output
461 : : //! surface) to side set ids read from the mesh file and removes those face
462 : : //! and node lists associated to side sets that the user did not set BCs or
463 : : //! listed as field output on (as they will not need processing further since
464 : : //! they will not be used).
465 : : //! \param[in,out] bnd Node or face lists mapped to side set ids
466 : : //! \return True if sidesets have been used and found in mesh
467 : : // *****************************************************************************
468 : : {
469 : : using PDETypes = ctr::parameters::Keys;
470 : : // Query side set ids at which BCs assigned for all BC types for all PDEs
471 : : using PDEsBCs =
472 : : tk::cartesian_product< PDETypes, ctr::bc::Keys >;
473 : 598 : std::unordered_set< int > usedsets;
474 [ + - ]: 299 : brigand::for_each< PDEsBCs >( UserBC( g_inputdeck, usedsets ) );
475 : :
476 : : // Query side sets of time dependent BCs (since tag::bctimedep is not a part
477 : : // of tag::bc)
478 [ + - ]: 299 : brigand::for_each< PDETypes >( UserTimedepBC(g_inputdeck, usedsets) );
479 : :
480 : : // Query side sets of boundaries prescribed as moving with ALE
481 [ + + ]: 301 : for (const auto& move : g_inputdeck.get< tag::ale, tag::move >())
482 [ + + ][ + - ]: 4 : for (auto i : move.get< tag::sideset >())
483 [ + - ][ + - ]: 2 : usedsets.insert( std::stoi(i) );
484 : :
485 : : // Add sidesets requested for field output
486 : 299 : const auto& ss = g_inputdeck.get< tag::cmd, tag::io, tag::surface >();
487 [ + + ]: 319 : for (const auto& s : ss) {
488 [ + - ]: 40 : std::stringstream conv( s );
489 : : int num;
490 [ + - ]: 20 : conv >> num;
491 [ + - ]: 20 : usedsets.insert( num );
492 : : }
493 : :
494 : : // Find user-configured side set ids among side sets read from mesh file
495 : 299 : std::unordered_set< int > sidesets_used;
496 [ + + ]: 1017 : for (auto i : usedsets) { // for all side sets used in control file
497 [ + - ][ + - ]: 718 : if (bnd.find(i) != end(bnd)) // used set found among side sets in file
498 [ + - ]: 718 : sidesets_used.insert( i ); // store side set id configured as BC
499 : : else {
500 [ - - ][ - - ]: 0 : Throw( "Boundary conditions specified on side set " +
[ - - ][ - - ]
[ - - ]
501 : : std::to_string(i) + " which does not exist in mesh file" );
502 : : }
503 : : }
504 : :
505 : : // Remove sidesets not used (will not process those further)
506 [ + - ]: 299 : tk::erase_if( bnd, [&]( auto& item ) {
507 [ + - ]: 828 : return sidesets_used.find( item.first ) == end(sidesets_used);
508 : : });
509 : :
510 : 598 : return !bnd.empty();
511 : : }
512 : :
513 : : void
514 : 233 : Transporter::createPartitioner()
515 : : // *****************************************************************************
516 : : // Create mesh partitioner AND boundary conditions group
517 : : // *****************************************************************************
518 : : {
519 : 233 : auto scheme = g_inputdeck.get< tag::discr, tag::scheme >();
520 [ + - ][ + - ]: 233 : auto centering = ctr::Scheme().centering( scheme );
521 [ + - ]: 466 : auto print = printer();
522 : :
523 : : // Create partitioner callbacks (order important)
524 : : tk::PartitionerCallback cbp {{
525 [ + - ][ + - ]: 233 : CkCallback( CkReductionTarget(Transporter,load), thisProxy )
526 [ + - ][ + - ]: 466 : , CkCallback( CkReductionTarget(Transporter,partitioned), thisProxy )
527 [ + - ][ + - ]: 466 : , CkCallback( CkReductionTarget(Transporter,distributed), thisProxy )
528 [ + - ][ + - ]: 466 : , CkCallback( CkReductionTarget(Transporter,refinserted), thisProxy )
529 [ + - ][ + - ]: 466 : , CkCallback( CkReductionTarget(Transporter,refined), thisProxy )
530 [ + - ]: 699 : }};
531 : :
532 : : // Create refiner callbacks (order important)
533 : : tk::RefinerCallback cbr {{
534 [ + - ][ + - ]: 233 : CkCallback( CkReductionTarget(Transporter,queriedRef), thisProxy )
535 [ + - ][ + - ]: 466 : , CkCallback( CkReductionTarget(Transporter,respondedRef), thisProxy )
536 [ + - ][ + - ]: 466 : , CkCallback( CkReductionTarget(Transporter,compatibility), thisProxy )
537 [ + - ][ + - ]: 466 : , CkCallback( CkReductionTarget(Transporter,bndint), thisProxy )
538 [ + - ][ + - ]: 466 : , CkCallback( CkReductionTarget(Transporter,matched), thisProxy )
539 [ + - ][ + - ]: 466 : , CkCallback( CkReductionTarget(Transporter,refined), thisProxy )
540 [ + - ]: 699 : }};
541 : :
542 : : // Create sorter callbacks (order important)
543 : : tk::SorterCallback cbs {{
544 [ + - ][ + - ]: 233 : CkCallback( CkReductionTarget(Transporter,queried), thisProxy )
545 [ + - ][ + - ]: 466 : , CkCallback( CkReductionTarget(Transporter,responded), thisProxy )
546 [ + - ][ + - ]: 466 : , CkCallback( CkReductionTarget(Transporter,discinserted), thisProxy )
547 [ + - ][ + - ]: 466 : , CkCallback( CkReductionTarget(Transporter,workinserted), thisProxy )
548 [ + - ]: 699 : }};
549 : :
550 : : // Start timer measuring preparation of mesh(es) for partitioning
551 [ + - ]: 233 : m_timer[ TimerTag::MESH_READ ];
552 : :
553 : : // Start preparing mesh(es)
554 [ + - ][ + - ]: 233 : print.diag( "Reading mesh(es)" );
555 : :
556 : : // Create (discretization) Scheme chare worker arrays for all meshes
557 [ + + ]: 466 : for ([[maybe_unused]] const auto& filename : m_input)
558 : : m_scheme.emplace_back( g_inputdeck.get< tag::discr, tag::scheme >(),
559 : : g_inputdeck.get< tag::ale, tag::ale >(),
560 [ + - ][ + - ]: 233 : need_linearsolver() );
561 : :
562 [ - + ][ - - ]: 233 : ErrChk( !m_input.empty(), "No input mesh" );
[ - - ][ - - ]
563 : :
564 : : // Read boundary (side set) data from a list of input mesh files
565 : 233 : std::size_t meshid = 0;
566 [ + + ]: 466 : for (const auto& filename : m_input) {
567 : : // Create mesh reader for reading side sets from file
568 [ + - ]: 466 : tk::MeshReader mr( filename );
569 : :
570 : : // Read out total number of mesh points from mesh file
571 [ + - ][ + - ]: 233 : m_npoin.push_back( mr.npoin() );
572 : :
573 : 466 : std::map< int, std::vector< std::size_t > > bface;
574 : 466 : std::map< int, std::vector< std::size_t > > faces;
575 : 466 : std::map< int, std::vector< std::size_t > > bnode;
576 : :
577 : : // Read boundary-face connectivity on side sets
578 [ + - ]: 233 : mr.readSidesetFaces( bface, faces );
579 : :
580 : 233 : bool bcs_set = false;
581 [ + + ]: 233 : if (centering == tk::Centering::ELEM) {
582 : :
583 : : // Verify boundarty condition (BC) side sets used exist in mesh file
584 [ + - ]: 96 : bcs_set = matchBCs( bface );
585 : :
586 [ + - ]: 137 : } else if (centering == tk::Centering::NODE) {
587 : :
588 : : // Read node lists on side sets
589 [ + - ]: 137 : bnode = mr.readSidesetNodes();
590 : : // Verify boundarty condition (BC) side sets used exist in mesh file
591 [ + - ]: 137 : bcs_set = matchBCs( bnode );
592 [ + + ][ + - ]: 137 : bcs_set = bcs_set || matchBCs( bface );
[ - + ]
593 : : }
594 : :
595 : : // Warn on no BCs
596 [ + + ][ + - ]: 233 : if (!bcs_set) print << "\n>>> WARNING: No boundary conditions set\n\n";
597 : :
598 [ + - ]: 233 : auto opt = m_scheme[meshid].arrayoptions();
599 : : // Create empty mesh refiner chare array (bound to workers)
600 [ + - ][ + - ]: 233 : m_refiner.push_back( CProxy_Refiner::ckNew(opt) );
[ + - ]
601 : : // Create empty mesh sorter Charm++ chare array (bound to workers)
602 [ + - ][ + - ]: 233 : m_sorter.push_back( CProxy_Sorter::ckNew(opt) );
[ + - ]
603 : :
604 : : // Create MeshWriter chare group for mesh
605 [ + - ][ + - ]: 233 : m_meshwriter.push_back(
606 : : tk::CProxy_MeshWriter::ckNew(
607 : 233 : g_inputdeck.get< tag::selected, tag::filetype >(),
608 : : centering,
609 : 233 : g_inputdeck.get< tag::cmd, tag::benchmark >(),
610 : 233 : m_input.size() ) );
611 : :
612 : : // Create mesh partitioner Charm++ chare nodegroup for all meshes
613 [ + - ]: 233 : m_partitioner.push_back(
614 : : CProxy_Partitioner::ckNew( meshid, filename, cbp, cbr, cbs,
615 [ + - ]: 233 : thisProxy, m_refiner.back(), m_sorter.back(), m_meshwriter.back(),
616 : 233 : m_scheme, bface, faces, bnode ) );
617 : :
618 : 233 : ++meshid;
619 : : }
620 : 233 : }
621 : :
622 : : void
623 : 233 : Transporter::load( std::size_t meshid, std::size_t nelem )
624 : : // *****************************************************************************
625 : : // Reduction target: the mesh has been read from file on all PEs
626 : : //! \param[in] meshid Mesh id (summed accross all compute nodes)
627 : : //! \param[in] nelem Number of mesh elements per mesh (summed across all
628 : : //! compute nodes)
629 : : // *****************************************************************************
630 : : {
631 : 233 : meshid /= static_cast< std::size_t >( CkNumNodes() );
632 [ - + ][ - - ]: 233 : Assert( meshid < m_nelem.size(), "MeshId indexing out" );
[ - - ][ - - ]
633 : 233 : m_nelem[meshid] = nelem;
634 : :
635 : : // Compute load distribution given total work (nelem) and user-specified
636 : : // virtualization
637 : : uint64_t chunksize, remainder;
638 : 466 : m_nchare[meshid] = static_cast<int>(
639 [ + - ]: 466 : tk::linearLoadDistributor(
640 : 233 : g_inputdeck.get< tag::cmd, tag::virtualization >(),
641 : 233 : m_nelem[meshid], CkNumPes(), chunksize, remainder ) );
642 : :
643 : : // Store sum of meshids (across all chares, key) for each meshid (value).
644 : : // This is used to look up the mesh id after collectives that sum their data.
645 [ + - ]: 233 : m_meshid[ static_cast<std::size_t>(m_nchare[meshid])*meshid ] = meshid;
646 [ - + ][ - - ]: 233 : Assert( meshid < m_nelem.size(), "MeshId indexing out" );
[ - - ][ - - ]
647 : :
648 : : // Partition first mesh
649 [ + - ][ + - ]: 233 : if (meshid == 0) m_partitioner[0].partition( m_nchare[0] );
650 : :
651 [ + - ]: 233 : if (++m_nload == m_nelem.size()) { // all meshes have been loaded
652 : 233 : m_nload = 0;
653 [ + - ]: 233 : auto print = printer();
654 : :
655 : : // Start timer measuring preparation of the mesh for partitioning
656 [ + - ]: 233 : const auto& timer = tk::cref_find( m_timer, TimerTag::MESH_READ );
657 [ + - ][ + - ]: 233 : print.diag( "Mesh read time: " + std::to_string( timer.dsec() ) + " sec" );
[ + - ][ + - ]
[ + - ]
658 : :
659 : : // Print out mesh partitioning configuration
660 [ + - ][ + - ]: 233 : print.section( "Mesh partitioning" );
661 : : print.Item< tk::ctr::PartitioningAlgorithm,
662 [ + - ]: 233 : tag::selected, tag::partitioner >();
663 [ + - ][ + - ]: 233 : print.item( "Virtualization [0.0...1.0]",
664 : 233 : g_inputdeck.get< tag::cmd, tag::virtualization >() );
665 : : // Print out initial mesh statistics
666 [ + - ][ + - ]: 233 : meshstat( "Initial load distribution" );
667 : :
668 : : // Tell meshwriter the total number of chares
669 [ + - ]: 233 : m_meshwriter[meshid].nchare( m_nchare[meshid] );
670 : :
671 : : // Query number of initial mesh refinement steps
672 : 233 : int nref = 0;
673 [ + + ]: 233 : if (g_inputdeck.get< tag::amr, tag::t0ref >())
674 : 25 : nref = static_cast<int>(g_inputdeck.get< tag::amr, tag::init >().size());
675 : :
676 [ + - ][ + - ]: 466 : m_progMesh.start( print, "Preparing mesh", {{ CkNumPes(), CkNumPes(), nref,
677 : 233 : m_nchare[0], m_nchare[0], m_nchare[0], m_nchare[0] }} );
678 : : }
679 : 233 : }
680 : :
681 : : void
682 : 233 : Transporter::partitioned( std::size_t meshid )
683 : : // *****************************************************************************
684 : : // Reduction target: a mesh has been partitioned
685 : : //! \param[in] meshid Mesh id
686 : : // *****************************************************************************
687 : : {
688 [ + - ]: 233 : if (++m_npart == m_nelem.size()) { // all meshes have been partitioned
689 : 233 : m_npart = 0;
690 : : } else // partition next mesh
691 : 0 : m_partitioner[meshid+1].partition( m_nchare[meshid+1] );
692 : 233 : }
693 : :
694 : : void
695 : 233 : Transporter::distributed( std::size_t meshid )
696 : : // *****************************************************************************
697 : : // Reduction target: all compute nodes have distributed their mesh after
698 : : // partitioning
699 : : //! \param[in] meshid Mesh id
700 : : // *****************************************************************************
701 : : {
702 : 233 : m_partitioner[meshid].refine();
703 : 233 : }
704 : :
705 : : void
706 : 233 : Transporter::refinserted( std::size_t meshid, std::size_t error )
707 : : // *****************************************************************************
708 : : // Reduction target: all compute nodes have created the mesh refiners
709 : : //! \param[in] meshid Mesh id (aggregated across all compute nodes with operator
710 : : //! max)
711 : : //! \param[in] error Error code (aggregated across all compute nodes with
712 : : //! operator max)
713 : : // *****************************************************************************
714 : : {
715 [ - + ]: 233 : if (error) {
716 : :
717 : 0 : printer() << "\n>>> ERROR: A worker chare was not assigned any mesh "
718 [ - - ][ - - ]: 0 : "elements after distributing mesh " + std::to_string(meshid) +
[ - - ]
719 : : ". This can happen in SMP-mode with a large +ppn "
720 : : "parameter (number of worker threads per logical node) and is "
721 : : "most likely the fault of the mesh partitioning algorithm not "
722 : : "tolerating the case when it is asked to divide the "
723 : : "computational domain into a number of partitions different "
724 : : "than the number of ranks it is called on, i.e., in case of "
725 : : "overdecomposition and/or calling the partitioner in SMP mode "
726 : : "with +ppn larger than 1. Solution 1: Try a different "
727 : : "partitioning algorithm (e.g., rcb instead of mj). Solution 2: "
728 [ - - ]: 0 : "Decrease +ppn.";
729 : 0 : finish( meshid );
730 : :
731 : : } else {
732 : :
733 : 233 : m_refiner[meshid].doneInserting();
734 : :
735 : : }
736 : 233 : }
737 : :
738 : : void
739 : 125 : Transporter::queriedRef( std::size_t meshid )
740 : : // *****************************************************************************
741 : : // Reduction target: all Refiner chares have queried their boundary edges
742 : : //! \param[in] meshid Mesh id
743 : : // *****************************************************************************
744 : : {
745 : 125 : m_refiner[meshid].response();
746 : 125 : }
747 : :
748 : : void
749 : 125 : Transporter::respondedRef( std::size_t meshid )
750 : : // *****************************************************************************
751 : : // Reduction target: all Refiner chares have setup their boundary edges
752 : : //! \param[in] meshid Mesh id
753 : : // *****************************************************************************
754 : : {
755 : 125 : m_refiner[meshid].refine();
756 : 125 : }
757 : :
758 : : void
759 : 78 : Transporter::compatibility( std::size_t meshid, std::size_t modified )
760 : : // *****************************************************************************
761 : : // Reduction target: all Refiner chares have received a round of edges,
762 : : // and have run their compatibility algorithm
763 : : //! \param[in] meshid Mesh id (aggregated across all chares using operator max)
764 : : //! \param[in] modified Modified flag, aggregated across all chares using
765 : : //! operator max), if nonzero, mesh is modified
766 : : //! \details This is called iteratively, until convergence by Refiner. At this
767 : : //! point all Refiner chares have received a round of edge data (tags whether
768 : : //! an edge needs to be refined, etc.), and applied the compatibility
769 : : //! algorithm independent of other Refiner chares. We keep going until the
770 : : //! mesh is no longer modified by the compatibility algorithm, based on a new
771 : : //! round of edge data communication started in Refiner::comExtra().
772 : : // *****************************************************************************
773 : : {
774 [ + + ]: 78 : if (modified)
775 : 2 : m_refiner[meshid].comExtra();
776 : : else
777 : 76 : m_refiner[meshid].correctref();
778 : 78 : }
779 : :
780 : : void
781 : 127 : Transporter::matched( std::size_t summeshid,
782 : : std::size_t nextra,
783 : : std::size_t nref,
784 : : std::size_t nderef,
785 : : std::size_t sumrefmode )
786 : : // *****************************************************************************
787 : : // Reduction target: all Refiner chares have matched/corrected the tagging
788 : : // of chare-boundary edges, all chares are ready to perform refinement.
789 : : //! \param[in] summeshid Mesh id (summed across all chares)
790 : : //! \param[in] nextra Sum (across all chares) of the number of edges on each
791 : : //! chare that need correction along chare boundaries
792 : : //! \param[in] nref Sum of number of refined tetrahedra across all chares.
793 : : //! \param[in] nderef Sum of number of derefined tetrahedra across all chares.
794 : : //! \param[in] sumrefmode Sum of contributions from all chares, encoding
795 : : //! refinement mode of operation.
796 : : // *****************************************************************************
797 : : {
798 : 127 : auto meshid = tk::cref_find( m_meshid, summeshid );
799 : :
800 : : // If at least a single edge on a chare still needs correction, do correction,
801 : : // otherwise, this mesh refinement step is complete
802 [ + + ]: 127 : if (nextra > 0) {
803 : :
804 : 2 : ++m_ncit[meshid];
805 : 2 : m_refiner[meshid].comExtra();
806 : :
807 : : } else {
808 : :
809 [ + - ]: 250 : auto print = printer();
810 : :
811 : : // decode refmode
812 : : auto refmode = static_cast< Refiner::RefMode >(
813 : 125 : sumrefmode / static_cast<std::size_t>(m_nchare[meshid]) );
814 : :
815 [ + + ]: 125 : if (refmode == Refiner::RefMode::T0REF) {
816 : :
817 [ + - ]: 66 : if (!g_inputdeck.get< tag::cmd, tag::feedback >()) {
818 [ + - ]: 66 : ctr::AMRInitial opt;
819 [ + - ][ + - ]: 990 : print.diag( { "meshid", "t0ref", "type", "nref", "nderef", "ncorr" },
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + + ]
[ + + ][ - - ]
[ - - ]
820 : : { std::to_string(meshid),
821 : 66 : std::to_string(m_nt0refit[meshid]),
822 : : "initial",
823 : : std::to_string(nref),
824 : : std::to_string(nderef),
825 : 858 : std::to_string(m_ncit[meshid]) } );
826 : 66 : ++m_nt0refit[meshid];
827 : : }
828 [ + - ]: 66 : m_progMesh.inc< REFINE >( print );
829 : :
830 [ + + ]: 59 : } else if (refmode == Refiner::RefMode::DTREF) {
831 : :
832 : 29 : auto dtref_uni = g_inputdeck.get< tag::amr, tag::dtref_uniform >();
833 [ + - ][ + - ]: 435 : print.diag( { "meshid", "dtref", "type", "nref", "nderef", "ncorr" },
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ][ + + ]
[ - - ][ - - ]
834 : : { std::to_string(meshid),
835 : 29 : std::to_string(++m_ndtrefit[meshid]),
836 : : (dtref_uni?"uniform":"error"),
837 : : std::to_string(nref),
838 : : std::to_string(nderef),
839 : 29 : std::to_string(m_ncit[meshid]) },
840 : 348 : false );
841 : :
842 [ + + ]: 30 : } else if (refmode == Refiner::RefMode::OUTREF) {
843 : :
844 [ + - ][ + - ]: 195 : print.diag( { "meshid", "outref", "nref", "nderef", "ncorr" },
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + + ]
[ + + ][ - - ]
[ - - ]
845 : : { std::to_string(meshid),
846 : 15 : std::to_string(++m_noutrefit[meshid]),
847 : : std::to_string(nref),
848 : : std::to_string(nderef),
849 : 165 : std::to_string(m_ncit[meshid]) }, false );
850 : :
851 [ + - ]: 15 : } else if (refmode == Refiner::RefMode::OUTDEREF) {
852 : :
853 [ + - ][ + - ]: 195 : print.diag( { "meshid", "outderef", "nref", "nderef", "ncorr" },
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + + ]
[ + + ][ - - ]
[ - - ]
854 : : { std::to_string(meshid),
855 : 15 : std::to_string(++m_noutderefit[meshid]),
856 : : std::to_string(nref),
857 : : std::to_string(nderef),
858 : 15 : std::to_string(m_ncit[meshid]) },
859 : 150 : false );
860 : :
861 [ - - ][ - - ]: 0 : } else Throw( "RefMode not implemented" );
[ - - ]
862 : :
863 : 125 : m_ncit[meshid] = 0;
864 [ + - ]: 125 : m_refiner[meshid].perform();
865 : :
866 : : }
867 : 127 : }
868 : :
869 : : void
870 : 221 : Transporter::bndint( tk::real sx, tk::real sy, tk::real sz, tk::real cb,
871 : : tk::real summeshid )
872 : : // *****************************************************************************
873 : : // Compute surface integral across the whole problem and perform leak-test
874 : : //! \param[in] sx X component of vector summed
875 : : //! \param[in] sy Y component of vector summed
876 : : //! \param[in] sz Z component of vector summed
877 : : //! \param[in] cb Invoke callback if positive
878 : : //! \param[in] summeshid Mesh id (summed accross all chares)
879 : : //! \details This function aggregates partial surface integrals across the
880 : : //! boundary faces of the whole problem. After this global sum a
881 : : //! non-zero vector result indicates a leak, e.g., a hole in the boundary,
882 : : //! which indicates an error in the boundary face data structures used to
883 : : //! compute the partial surface integrals.
884 : : // *****************************************************************************
885 : : {
886 [ + - ]: 221 : auto meshid = tk::cref_find( m_meshid, static_cast<std::size_t>(summeshid) );
887 : :
888 [ + - ]: 442 : std::stringstream err;
889 [ + + ]: 221 : if (cb < 0.0) {
890 : : err << "Mesh boundary leaky after mesh refinement step; this is due to a "
891 : : "problem with updating the side sets used to specify boundary conditions "
892 [ + - ]: 125 : "on faces: ";
893 [ + - ]: 96 : } else if (cb > 0.0) {
894 : : err << "Mesh boundary leaky during initialization; this is due to "
895 : : "incorrect or incompletely specified boundary conditions for a given input "
896 [ + - ]: 96 : "mesh: ";
897 : : }
898 : :
899 : 221 : auto eps = 1.0e-10;
900 [ + - ][ + - ]: 221 : if (std::abs(sx) > eps || std::abs(sy) > eps || std::abs(sz) > eps) {
[ - + ][ - + ]
901 [ - - ][ - - ]: 0 : err << "Integral result must be a zero vector: " << std::setprecision(12) <<
902 [ - - ][ - - ]: 0 : std::abs(sx) << ", " << std::abs(sy) << ", " << std::abs(sz) <<
[ - - ][ - - ]
[ - - ]
903 [ - - ][ - - ]: 0 : ", eps = " << eps;
904 [ - - ][ - - ]: 0 : Throw( err.str() );
[ - - ]
905 : : }
906 : :
907 [ + + ][ + - ]: 221 : if (cb > 0.0) m_scheme[meshid].bcast< Scheme::resizeComm >();
908 : 221 : }
909 : :
910 : : void
911 : 233 : Transporter::refined( std::size_t summeshid,
912 : : std::size_t nelem,
913 : : std::size_t npoin )
914 : : // *****************************************************************************
915 : : // Reduction target: all chares have refined their mesh
916 : : //! \param[in] summeshid Mesh id (summed accross all Refiner chares)
917 : : //! \param[in] nelem Total number of elements in mesh summed across the
918 : : //! distributed mesh
919 : : //! \param[in] npoin Total number of mesh points summed across the distributed
920 : : //! mesh. Note that in parallel this is larger than the number of points in
921 : : //! the mesh, because the boundary nodes are multi-counted. But we only need
922 : : //! an equal or larger than npoin for Sorter::setup, so this is okay.
923 : : // *****************************************************************************
924 : : {
925 : 233 : auto meshid = tk::cref_find( m_meshid, summeshid );
926 : :
927 : : // Store new number of elements for initially refined mesh
928 : 233 : m_nelem[meshid] = nelem;
929 : :
930 : 233 : m_sorter[meshid].doneInserting();
931 : 233 : m_sorter[meshid].setup( npoin );
932 : 233 : }
933 : :
934 : : void
935 : 233 : Transporter::queried( std::size_t meshid )
936 : : // *****************************************************************************
937 : : // Reduction target: all Sorter chares have queried their boundary edges
938 : : //! \param[in] meshid Mesh id
939 : : // *****************************************************************************
940 : : {
941 : 233 : m_sorter[meshid].response();
942 : 233 : }
943 : :
944 : : void
945 : 233 : Transporter::responded( std::size_t meshid )
946 : : // *****************************************************************************
947 : : // Reduction target: all Sorter chares have responded with their boundary edges
948 : : //! \param[in] meshid Mesh id
949 : : // *****************************************************************************
950 : : {
951 : 233 : m_sorter[meshid].start();
952 : 233 : }
953 : :
954 : : void
955 : 373 : Transporter::resized( std::size_t meshid )
956 : : // *****************************************************************************
957 : : // Reduction target: all worker chares have resized their own mesh data after
958 : : // AMR or ALE
959 : : //! \param[in] meshid Mesh id
960 : : //! \note Only used for nodal schemes
961 : : // *****************************************************************************
962 : : {
963 : 373 : m_scheme[meshid].disc().vol();
964 : 373 : m_scheme[meshid].bcast< Scheme::lhs >();
965 : 373 : }
966 : :
967 : : void
968 : 112 : Transporter::startEsup( std::size_t meshid )
969 : : // *****************************************************************************
970 : : // Reduction target: all worker chares have generated their own esup
971 : : //! \param[in] meshid Mesh id
972 : : //! \note Only used for cell-centered schemes
973 : : // *****************************************************************************
974 : : {
975 : 112 : m_scheme[meshid].bcast< Scheme::nodeNeighSetup >();
976 : 112 : }
977 : :
978 : : void
979 : 233 : Transporter::discinserted( std::size_t meshid )
980 : : // *****************************************************************************
981 : : // Reduction target: all Discretization chares have been inserted
982 : : //! \param[in] meshid Mesh id
983 : : // *****************************************************************************
984 : : {
985 : 233 : m_scheme[meshid].disc().doneInserting();
986 : 233 : }
987 : :
988 : : void
989 : 258 : Transporter::meshstat( const std::string& header ) const
990 : : // *****************************************************************************
991 : : // Print out mesh statistics
992 : : //! \param[in] header Section header
993 : : // *****************************************************************************
994 : : {
995 [ + - ]: 516 : auto print = printer();
996 : :
997 [ + - ]: 258 : print.section( header );
998 : :
999 [ - + ]: 258 : if (m_nelem.size() > 1) {
1000 [ - - ][ - - ]: 0 : print.item( "Number of tetrahedra (per mesh)",tk::parameters(m_nelem) );
[ - - ]
1001 [ - - ][ - - ]: 0 : print.item( "Number of points (per mesh)", tk::parameters(m_npoin) );
[ - - ]
1002 [ - - ][ - - ]: 0 : print.item( "Number of work units (per mesh)", tk::parameters(m_nchare) );
[ - - ]
1003 : : }
1004 : :
1005 [ + - ][ + - ]: 258 : print.item( "Total number of tetrahedra",
1006 : 258 : std::accumulate( begin(m_nelem), end(m_nelem), 0UL ) );
1007 [ + - ][ + - ]: 258 : print.item( "Total number of points",
1008 : 258 : std::accumulate( begin(m_npoin), end(m_npoin), 0UL ) );
1009 [ + - ][ + - ]: 258 : print.item( "Total number of work units",
1010 : 258 : std::accumulate( begin(m_nchare), end(m_nchare), 0 ) );
1011 : :
1012 [ + - ]: 258 : print.endsubsection();
1013 : 258 : }
1014 : :
1015 : : bool
1016 : 466 : Transporter::need_linearsolver() const
1017 : : // *****************************************************************************
1018 : : // Decide if we need a linear solver for ALE
1019 : : //! \return True if ALE will neeed a linear solver
1020 : : // *****************************************************************************
1021 : : {
1022 : 466 : auto smoother = g_inputdeck.get< tag::ale, tag::smoother >();
1023 : :
1024 [ + + ][ + + ]: 476 : if ( g_inputdeck.get< tag::ale, tag::ale >() and
[ + + ]
1025 [ + + ]: 10 : (smoother == ctr::MeshVelocitySmootherType::LAPLACE ||
1026 : : smoother == ctr::MeshVelocitySmootherType::HELMHOLTZ) )
1027 : : {
1028 : 12 : return true;
1029 : : } else {
1030 : 454 : return false;
1031 : : }
1032 : : }
1033 : :
1034 : : void
1035 : 233 : Transporter::disccreated( std::size_t summeshid, std::size_t npoin )
1036 : : // *****************************************************************************
1037 : : // Reduction target: all Discretization constructors have been called
1038 : : //! \param[in] summeshid Mesh id (summed accross all chares)
1039 : : //! \param[in] npoin Total number of mesh points (summed across all chares)
1040 : : //! Note that as opposed to npoin in refined(), this npoin is not
1041 : : //! multi-counted, and thus should be correct in parallel.
1042 : : // *****************************************************************************
1043 : : {
1044 : 233 : auto meshid = tk::cref_find( m_meshid, summeshid );
1045 : :
1046 : : // Update number of mesh points for mesh, since it may have been refined
1047 [ + + ]: 233 : if (g_inputdeck.get< tag::amr, tag::t0ref >()) m_npoin[meshid] = npoin;
1048 : :
1049 [ + - ]: 233 : if (++m_ndisc == m_nelem.size()) { // all Disc arrays have been created
1050 : 233 : m_ndisc = 0;
1051 [ + - ]: 466 : auto print = printer();
1052 [ + - ]: 233 : m_progMesh.end( print );
1053 [ + + ]: 233 : if (g_inputdeck.get< tag::amr, tag::t0ref >())
1054 [ + - ][ + - ]: 25 : meshstat( "Initially (t<0) refined mesh graph statistics" );
1055 : : }
1056 : :
1057 : 233 : m_refiner[meshid].sendProxy();
1058 : :
1059 [ + + ]: 233 : if (g_inputdeck.get< tag::discr, tag::scheme >() == ctr::SchemeType::DiagCG)
1060 : 83 : m_scheme[meshid].fct().doneInserting();
1061 : :
1062 [ + + ]: 233 : if (g_inputdeck.get< tag::ale, tag::ale >())
1063 : 10 : m_scheme[meshid].ale().doneInserting();
1064 : :
1065 [ + + ]: 233 : if (need_linearsolver())
1066 : 6 : m_scheme[meshid].conjugategradients().doneInserting();
1067 : :
1068 : 233 : m_scheme[meshid].disc().vol();
1069 : 233 : }
1070 : :
1071 : : void
1072 : 233 : Transporter::workinserted( std::size_t meshid )
1073 : : // *****************************************************************************
1074 : : // Reduction target: all worker (derived discretization) chares have been
1075 : : // inserted
1076 : : //! \param[in] meshid Mesh id
1077 : : // *****************************************************************************
1078 : : {
1079 : 233 : m_scheme[meshid].bcast< Scheme::doneInserting >();
1080 : 233 : }
1081 : :
1082 : : void
1083 : 233 : Transporter::diagHeader()
1084 : : // *****************************************************************************
1085 : : // Configure and write diagnostics file header
1086 : : // *****************************************************************************
1087 : : {
1088 : : // Output header for diagnostics output file
1089 : 233 : tk::DiagWriter dw( g_inputdeck.get< tag::cmd, tag::io, tag::diag >(),
1090 : 233 : g_inputdeck.get< tag::flformat, tag::diag >(),
1091 [ + - ]: 466 : g_inputdeck.get< tag::prec, tag::diag >() );
1092 : :
1093 : : // Collect variables names for integral/diagnostics output
1094 : 466 : std::vector< std::string > var;
1095 : 233 : const auto scheme = g_inputdeck.get< tag::discr, tag::scheme >();
1096 [ + + ][ + + ]: 233 : if (scheme == ctr::SchemeType::DiagCG || scheme == ctr::SchemeType::ALECG)
1097 [ + + ][ + - ]: 274 : for (const auto& eq : g_cgpde) varnames( eq, var );
1098 [ + + ][ + + ]: 96 : else if (scheme == ctr::SchemeType::DG ||
1099 [ + + ][ + + ]: 38 : scheme == ctr::SchemeType::P0P1 || scheme == ctr::SchemeType::DGP1 ||
1100 [ + - ]: 12 : scheme == ctr::SchemeType::DGP2 || scheme == ctr::SchemeType::PDG)
1101 [ + + ][ + - ]: 192 : for (const auto& eq : g_dgpde) varnames( eq, var );
1102 [ - - ][ - - ]: 0 : else Throw( "Diagnostics header not handled for discretization scheme" );
[ - - ]
1103 : :
1104 [ + - ]: 466 : const tk::ctr::Error opt;
1105 : 233 : auto nv = var.size();
1106 : 466 : std::vector< std::string > d;
1107 : :
1108 : : // Add 'L2(var)' for all variables as those are always computed
1109 [ + - ]: 233 : const auto& l2name = opt.name( tk::ctr::ErrorType::L2 );
1110 [ + + ][ + - ]: 895 : for (std::size_t i=0; i<nv; ++i) d.push_back( l2name + '(' + var[i] + ')' );
[ + - ][ + - ]
[ + - ]
1111 : :
1112 : : // Query user-requested diagnostics and augment diagnostics file header by
1113 : : // 'err(var)', where 'err' is the error type configured, and var is the
1114 : : // variable computed, for all variables and all error types configured.
1115 : 233 : const auto& err = g_inputdeck.get< tag::diag, tag::error >();
1116 [ + + ]: 381 : for (const auto& e : err) {
1117 [ + - ]: 148 : const auto& errname = opt.name( e );
1118 [ + + ]: 730 : for (std::size_t i=0; i<nv; ++i)
1119 [ + - ][ + - ]: 582 : d.push_back( errname + '(' + var[i] + "-IC)" );
[ + - ][ + - ]
1120 : : }
1121 : :
1122 : : // Augment diagnostics variables by L2-norm of the residual and total energy
1123 [ + + ][ + + ]: 233 : if (scheme == ctr::SchemeType::DiagCG || scheme == ctr::SchemeType::ALECG) {
1124 [ + + ][ + - ]: 515 : for (std::size_t i=0; i<nv; ++i) d.push_back( "L2(d" + var[i] + ')' );
[ + - ][ + - ]
1125 : : }
1126 [ + - ][ + - ]: 233 : d.push_back( "mE" );
1127 : :
1128 : : // Write diagnostics header
1129 [ + - ]: 233 : dw.header( d );
1130 : 233 : }
1131 : :
1132 : : void
1133 : 249 : Transporter::comfinal( std::size_t initial, std::size_t summeshid )
1134 : : // *****************************************************************************
1135 : : // Reduction target indicating that communication maps have been setup
1136 : : //! \param[in] initial Sum of contributions from all chares. If larger than
1137 : : //! zero, we are during time stepping and if zero we are during setup.
1138 : : //! \param[in] summeshid Mesh id (summed accross the distributed mesh)
1139 : : // *****************************************************************************
1140 : : // [Discretization-specific communication maps]
1141 : : {
1142 [ + - ]: 249 : auto meshid = tk::cref_find( m_meshid, static_cast<std::size_t>(summeshid) );
1143 : :
1144 [ + + ]: 249 : if (initial > 0) {
1145 : 233 : m_scheme[meshid].bcast< Scheme::setup >();
1146 : : // Turn on automatic load balancing
1147 [ + - ]: 233 : if (++m_ncom == m_nelem.size()) { // all worker arrays have finished
1148 : 233 : m_ncom = 0;
1149 [ + - ]: 233 : auto print = printer();
1150 [ + - ]: 233 : m_progWork.end( print );
1151 [ + - ]: 233 : tk::CProxy_LBSwitch::ckNew();
1152 [ + - ][ + - ]: 233 : print.diag( "Load balancing on (if enabled in Charm++)" );
1153 : : }
1154 : : } else {
1155 : 16 : m_scheme[meshid].bcast< Scheme::lhs >();
1156 : : }
1157 : 249 : }
1158 : : // [Discretization-specific communication maps]
1159 : :
1160 : : void
1161 : 606 : Transporter::totalvol( tk::real v, tk::real initial, tk::real summeshid )
1162 : : // *****************************************************************************
1163 : : // Reduction target summing total mesh volume across all workers
1164 : : //! \param[in] v Mesh volume summed across the distributed mesh
1165 : : //! \param[in] initial Sum of contributions from all chares. If larger than
1166 : : //! zero, we are during setup, if zero, during time stepping.
1167 : : //! \param[in] summeshid Mesh id (summed accross the distributed mesh)
1168 : : // *****************************************************************************
1169 : : {
1170 [ + - ]: 606 : auto meshid = tk::cref_find( m_meshid, static_cast<std::size_t>(summeshid) );
1171 : :
1172 : 606 : m_meshvol[meshid] = v;
1173 : :
1174 [ + + ]: 606 : if (initial > 0.0) // during initialization
1175 : 233 : m_scheme[meshid].disc().stat( v );
1176 : : else // during ALE or AMR
1177 : 373 : m_scheme[meshid].bcast< Scheme::resized >();
1178 : 606 : }
1179 : :
1180 : : void
1181 : 233 : Transporter::minstat( tk::real d0, tk::real d1, tk::real d2, tk::real rmeshid )
1182 : : // *****************************************************************************
1183 : : // Reduction target yielding minimum mesh statistcs across all workers
1184 : : //! \param[in] d0 Minimum mesh statistics collected over all chares
1185 : : //! \param[in] d1 Minimum mesh statistics collected over all chares
1186 : : //! \param[in] d2 Minimum mesh statistics collected over all chares
1187 : : //! \param[in] rmeshid Mesh id as a real
1188 : : // *****************************************************************************
1189 : : {
1190 : 233 : auto meshid = static_cast<std::size_t>(rmeshid);
1191 : :
1192 : 233 : m_minstat[meshid][0] = d0; // minimum edge length
1193 : 233 : m_minstat[meshid][1] = d1; // minimum cell volume cubic root
1194 : 233 : m_minstat[meshid][2] = d2; // minimum number of cells on chare
1195 : :
1196 : 233 : minstat_complete(meshid);
1197 : 233 : }
1198 : :
1199 : : void
1200 : 233 : Transporter::maxstat( tk::real d0, tk::real d1, tk::real d2, tk::real rmeshid )
1201 : : // *****************************************************************************
1202 : : // Reduction target yielding the maximum mesh statistics across all workers
1203 : : //! \param[in] d0 Maximum mesh statistics collected over all chares
1204 : : //! \param[in] d1 Maximum mesh statistics collected over all chares
1205 : : //! \param[in] d2 Maximum mesh statistics collected over all chares
1206 : : //! \param[in] rmeshid Mesh id as a real
1207 : : // *****************************************************************************
1208 : : {
1209 : 233 : auto meshid = static_cast<std::size_t>(rmeshid);
1210 : :
1211 : 233 : m_maxstat[meshid][0] = d0; // maximum edge length
1212 : 233 : m_maxstat[meshid][1] = d1; // maximum cell volume cubic root
1213 : 233 : m_maxstat[meshid][2] = d2; // maximum number of cells on chare
1214 : :
1215 : 233 : maxstat_complete(meshid);
1216 : 233 : }
1217 : :
1218 : : void
1219 : 233 : Transporter::sumstat( tk::real d0, tk::real d1, tk::real d2, tk::real d3,
1220 : : tk::real d4, tk::real d5, tk::real summeshid )
1221 : : // *****************************************************************************
1222 : : // Reduction target yielding the sum mesh statistics across all workers
1223 : : //! \param[in] d0 Sum mesh statistics collected over all chares
1224 : : //! \param[in] d1 Sum mesh statistics collected over all chares
1225 : : //! \param[in] d2 Sum mesh statistics collected over all chares
1226 : : //! \param[in] d3 Sum mesh statistics collected over all chares
1227 : : //! \param[in] d4 Sum mesh statistics collected over all chares
1228 : : //! \param[in] d5 Sum mesh statistics collected over all chares
1229 : : //! \param[in] summeshid Mesh id (summed accross the distributed mesh)
1230 : : // *****************************************************************************
1231 : : {
1232 [ + - ]: 233 : auto meshid = tk::cref_find( m_meshid, static_cast<std::size_t>(summeshid) );
1233 : :
1234 : 233 : m_avgstat[meshid][0] = d1 / d0; // average edge length
1235 : 233 : m_avgstat[meshid][1] = d3 / d2; // average cell volume cubic root
1236 : 233 : m_avgstat[meshid][2] = d5 / d4; // average number of cells per chare
1237 : :
1238 : 233 : sumstat_complete(meshid);
1239 : 233 : }
1240 : :
1241 : : void
1242 : 233 : Transporter::pdfstat( CkReductionMsg* msg )
1243 : : // *****************************************************************************
1244 : : // Reduction target yielding PDF of mesh statistics across all workers
1245 : : //! \param[in] msg Serialized PDF
1246 : : // *****************************************************************************
1247 : : {
1248 : : std::size_t meshid;
1249 : 466 : std::vector< tk::UniPDF > pdf;
1250 : :
1251 : : // Deserialize final PDF
1252 : 466 : PUP::fromMem creator( msg->getData() );
1253 [ + - ]: 233 : creator | meshid;
1254 [ + - ]: 233 : creator | pdf;
1255 [ + - ][ + - ]: 233 : delete msg;
1256 : :
1257 [ + - ]: 466 : auto id = std::to_string(meshid);
1258 : :
1259 : : // Create new PDF file (overwrite if exists)
1260 [ + - ][ + - ]: 699 : tk::PDFWriter pdfe( "mesh_edge_pdf." + id + ".txt" );
[ + - ]
1261 : : // Output edgelength PDF
1262 : : // cppcheck-suppress containerOutOfBounds
1263 [ + - ]: 233 : pdfe.writeTxt( pdf[0],
1264 [ + - ][ + - ]: 699 : tk::ctr::PDFInfo{ {"PDF"}, {}, {"edgelength"}, 0, 0.0 } );
[ + - ][ + + ]
[ - - ]
1265 : :
1266 : : // Create new PDF file (overwrite if exists)
1267 [ + - ][ + - ]: 699 : tk::PDFWriter pdfv( "mesh_vol_pdf." + id + ".txt" );
[ + - ]
1268 : : // Output cell volume cubic root PDF
1269 [ + - ]: 233 : pdfv.writeTxt( pdf[1],
1270 [ + - ][ + - ]: 699 : tk::ctr::PDFInfo{ {"PDF"}, {}, {"V^{1/3}"}, 0, 0.0 } );
[ + - ][ + + ]
[ - - ]
1271 : :
1272 : : // Create new PDF file (overwrite if exists)
1273 [ + - ][ + - ]: 699 : tk::PDFWriter pdfn( "mesh_ntet_pdf." + id + ".txt" );
[ + - ]
1274 : : // Output number of cells PDF
1275 [ + - ]: 233 : pdfn.writeTxt( pdf[2],
1276 [ + - ][ + - ]: 699 : tk::ctr::PDFInfo{ {"PDF"}, {}, {"ntets"}, 0, 0.0 } );
[ + - ][ + + ]
[ - - ]
1277 : :
1278 [ + - ]: 233 : pdfstat_complete(meshid);
1279 : 233 : }
1280 : :
1281 : : void
1282 : 233 : Transporter::stat()
1283 : : // *****************************************************************************
1284 : : // Echo diagnostics on mesh statistics
1285 : : // *****************************************************************************
1286 : : {
1287 [ + - ]: 466 : auto print = printer();
1288 : :
1289 [ + - ]: 233 : if (++m_nstat == m_nelem.size()) { // stats from all meshes have arrived
1290 : 233 : m_nstat = 0;
1291 [ + + ]: 466 : for (std::size_t i=0; i<m_nelem.size(); ++i) {
1292 [ + - ]: 233 : print.diag(
1293 [ + - ][ + - ]: 466 : "Mesh " + std::to_string(i) +
[ + - ]
1294 [ + - ]: 466 : " distribution statistics: min/max/avg(edgelength) = " +
1295 [ + - ][ + - ]: 932 : std::to_string( m_minstat[i][0] ) + " / " +
[ + - ]
1296 [ + - ][ + - ]: 932 : std::to_string( m_maxstat[i][0] ) + " / " +
[ + - ]
1297 [ + - ][ + - ]: 932 : std::to_string( m_avgstat[i][0] ) + ", " +
[ + - ]
1298 [ + - ]: 466 : "min/max/avg(V^{1/3}) = " +
1299 [ + - ][ + - ]: 932 : std::to_string( m_minstat[i][1] ) + " / " +
[ + - ]
1300 [ + - ][ + - ]: 932 : std::to_string( m_maxstat[i][1] ) + " / " +
[ + - ]
1301 [ + - ][ + - ]: 932 : std::to_string( m_avgstat[i][1] ) + ", " +
[ + - ]
1302 [ + - ]: 466 : "min/max/avg(ntets) = " +
1303 [ + - ][ + - ]: 932 : std::to_string( static_cast<std::size_t>(m_minstat[i][2]) ) + " / " +
[ + - ]
1304 [ + - ][ + - ]: 932 : std::to_string( static_cast<std::size_t>(m_maxstat[i][2]) ) + " / " +
[ + - ]
1305 [ + - ]: 466 : std::to_string( static_cast<std::size_t>(m_avgstat[i][2]) ) );
1306 : : }
1307 : :
1308 : : // Print out time integration header to screen
1309 [ + - ]: 233 : inthead( print );
1310 : :
1311 [ + - ][ + - ]: 233 : m_progWork.start( print, "Preparing workers",
1312 : 233 : {{ m_nchare[0], m_nchare[0], m_nchare[0], m_nchare[0], m_nchare[0] }} );
1313 : :
1314 : : // Create "derived-class" workers
1315 [ + + ][ + - ]: 466 : for (std::size_t i=0; i<m_nelem.size(); ++i) m_sorter[i].createWorkers();
1316 : : }
1317 : 233 : }
1318 : :
1319 : : void
1320 : 233 : Transporter::boxvol( tk::real v, tk::real summeshid )
1321 : : // *****************************************************************************
1322 : : // Reduction target computing total volume of IC box
1323 : : //! \param[in] v Total volume within user-specified box IC
1324 : : //! \param[in] summeshid Mesh id as a real (summed accross the distributed mesh)
1325 : : // *****************************************************************************
1326 : : {
1327 [ + - ]: 233 : auto meshid = tk::cref_find( m_meshid, static_cast<std::size_t>(summeshid) );
1328 [ + + ][ + - ]: 233 : if (v > 0.0) printer().diag( "Box IC volume: " + std::to_string(v) );
[ + - ][ + - ]
1329 : 233 : m_scheme[meshid].bcast< Scheme::box >( v );
1330 : 233 : }
1331 : :
1332 : : void
1333 : 233 : Transporter::inthead( const InciterPrint& print )
1334 : : // *****************************************************************************
1335 : : // Print out time integration header to screen
1336 : : //! \param[in] print Pretty printer object to use for printing
1337 : : // *****************************************************************************
1338 : : {
1339 : 233 : auto refined = g_inputdeck.get< tag::cmd, tag::io, tag::refined >();
1340 : 233 : const auto scheme = g_inputdeck.get< tag::discr, tag::scheme >();
1341 [ + + ][ - + ]: 233 : if (refined && scheme == ctr::SchemeType::DG) {
1342 [ - - ]: 0 : printer() << "\n>>> WARNING: Ignoring refined field output for DG(P0)\n\n";
1343 : 0 : refined = false;
1344 : : }
1345 : :
1346 [ + - ][ + - ]: 699 : print.inthead( "Time integration", "Navier-Stokes solver",
[ + - ][ + - ]
1347 : : "Legend: it - iteration count\n"
1348 : : " t - physics time\n"
1349 : : " dt - physics time step size\n"
1350 : : " ETE - estimated wall-clock time elapsed (h:m:s)\n"
1351 : : " ETA - estimated wall-clock time for accomplishment (h:m:s)\n"
1352 : : " EGT - estimated grind wall-clock time (ms/timestep)\n"
1353 : : " flg - status flags, legend:\n"
1354 [ + + ][ + - ]: 466 : " f - " + std::string(refined ? "refined " : "")
[ + - ]
1355 [ + - ]: 466 : + "field (volume and surface)\n"
1356 : : " d - diagnostics\n"
1357 : : " t - physics time history\n"
1358 : : " h - h-refinement\n"
1359 : : " l - load balancing\n"
1360 : : " r - checkpoint\n"
1361 : : " a - ALE mesh velocity linear solver did not converge\n",
1362 : : "\n it t dt ETE ETA EGT flg\n"
1363 : : " -------------------------------------------------------------------------\n" );
1364 : 233 : }
1365 : :
1366 : : void
1367 : 5041 : Transporter::diagnostics( CkReductionMsg* msg )
1368 : : // *****************************************************************************
1369 : : // Reduction target optionally collecting diagnostics, e.g., residuals
1370 : : //! \param[in] msg Serialized diagnostics vector aggregated across all PEs
1371 : : //! \note Only used for nodal schemes
1372 : : // *****************************************************************************
1373 : : {
1374 : : std::size_t meshid;
1375 : 10082 : std::vector< std::vector< tk::real > > d;
1376 : :
1377 : : // Deserialize diagnostics vector
1378 : 10082 : PUP::fromMem creator( msg->getData() );
1379 [ + - ]: 5041 : creator | meshid;
1380 [ + - ]: 5041 : creator | d;
1381 [ + - ][ + - ]: 5041 : delete msg;
1382 : :
1383 [ + - ]: 10082 : auto id = std::to_string(meshid);
1384 : 5041 : auto ncomp = g_inputdeck.get< tag::component >().nprop();
1385 : :
1386 [ - + ][ - - ]: 5041 : Assert( ncomp > 0, "Number of scalar components must be positive");
[ - - ][ - - ]
1387 [ - + ][ - - ]: 5041 : Assert( d.size() == NUMDIAG, "Diagnostics vector size mismatch" );
[ - - ][ - - ]
1388 : :
1389 [ + + ]: 45369 : for (std::size_t i=0; i<d.size(); ++i)
1390 [ - + ][ - - ]: 40328 : Assert( d[i].size() == ncomp, "Size mismatch at final stage of "
[ - - ][ - - ]
1391 : : "diagnostics aggregation for mesh " + id );
1392 : :
1393 : : // Allocate storage for those diagnostics that are always computed
1394 [ + - ]: 10082 : std::vector< tk::real > diag( ncomp, 0.0 );
1395 : :
1396 : : // Finish computing diagnostics
1397 [ + + ]: 26617 : for (std::size_t i=0; i<d[L2SOL].size(); ++i)
1398 : 21576 : diag[i] = sqrt( d[L2SOL][i] / m_meshvol[meshid] );
1399 : :
1400 : : // Query user-requested error types to output
1401 : 5041 : const auto& error = g_inputdeck.get< tag::diag, tag::error >();
1402 : :
1403 : 5041 : decltype(ncomp) n = 0;
1404 [ + + ]: 9290 : for (const auto& e : error) {
1405 : 4249 : n += ncomp;
1406 [ + + ]: 4249 : if (e == tk::ctr::ErrorType::L2) {
1407 : : // Finish computing the L2 norm of the numerical - analytical solution
1408 [ + + ]: 24421 : for (std::size_t i=0; i<d[L2ERR].size(); ++i)
1409 [ + - ]: 20300 : diag.push_back( sqrt( d[L2ERR][i] / m_meshvol[meshid] ) );
1410 [ + - ]: 128 : } else if (e == tk::ctr::ErrorType::LINF) {
1411 : : // Finish computing the Linf norm of the numerical - analytical solution
1412 [ + + ]: 659 : for (std::size_t i=0; i<d[LINFERR].size(); ++i)
1413 [ + - ]: 531 : diag.push_back( d[LINFERR][i] );
1414 : : }
1415 : : }
1416 : :
1417 : : // Finish computing the L2 norm of the residual and append
1418 : 5041 : const auto scheme = g_inputdeck.get< tag::discr, tag::scheme >();
1419 [ + - ]: 10082 : std::vector< tk::real > l2res( d[L2RES].size(), 0.0 );
1420 [ + + ][ + + ]: 5041 : if (scheme == ctr::SchemeType::DiagCG || scheme == ctr::SchemeType::ALECG) {
1421 [ + + ]: 16043 : for (std::size_t i=0; i<d[L2RES].size(); ++i) {
1422 : 12807 : l2res[i] = std::sqrt( d[L2RES][i] / m_meshvol[meshid] );
1423 [ + - ]: 12807 : diag.push_back( l2res[i] );
1424 : : }
1425 : : }
1426 : :
1427 : : // Append total energy
1428 [ + - ]: 5041 : diag.push_back( d[TOTALSOL][0] );
1429 : :
1430 : : // Append diagnostics file at selected times
1431 [ + - ]: 10082 : auto filename = g_inputdeck.get< tag::cmd, tag::io, tag::diag >();
1432 [ - + ][ - - ]: 5041 : if (m_nelem.size() > 1) filename += '.' + id;
[ - - ]
1433 : : tk::DiagWriter dw( filename,
1434 : 5041 : g_inputdeck.get< tag::flformat, tag::diag >(),
1435 : 5041 : g_inputdeck.get< tag::prec, tag::diag >(),
1436 [ + - ]: 10082 : std::ios_base::app );
1437 [ + - ]: 5041 : dw.diag( static_cast<uint64_t>(d[ITER][0]), d[TIME][0], d[DT][0], diag );
1438 : :
1439 : : // Continue time step
1440 [ + - ]: 5041 : m_scheme[meshid].bcast< Scheme::refine >( l2res );
1441 : 5041 : }
1442 : :
1443 : : void
1444 : 233 : Transporter::resume()
1445 : : // *****************************************************************************
1446 : : // Resume execution from checkpoint/restart files
1447 : : //! \details This is invoked by Charm++ after the checkpoint is done, as well as
1448 : : //! when the restart (returning from a checkpoint) is complete
1449 : : // *****************************************************************************
1450 : : {
1451 [ - + ]: 466 : if (std::any_of(begin(m_finished), end(m_finished), [](auto f){return !f;})) {
1452 : : // If just restarted from a checkpoint, Main( CkMigrateMessage* msg ) has
1453 : : // increased nrestart in g_inputdeck, but only on PE 0, so broadcast.
1454 : 0 : auto nrestart = g_inputdeck.get< tag::cmd, tag::io, tag::nrestart >();
1455 [ - - ]: 0 : for (std::size_t i=0; i<m_nelem.size(); ++i)
1456 [ - - ]: 0 : m_scheme[i].bcast< Scheme::evalLB >( nrestart );
1457 : : } else
1458 : 233 : mainProxy.finalize();
1459 : 233 : }
1460 : :
1461 : : void
1462 : 233 : Transporter::checkpoint( std::size_t finished, std::size_t meshid )
1463 : : // *****************************************************************************
1464 : : // Save checkpoint/restart files
1465 : : //! \param[in] finished Nonzero if finished with time stepping
1466 : : //! \param[in] meshid Mesh id
1467 : : // *****************************************************************************
1468 : : {
1469 : 233 : m_finished[meshid] = finished;
1470 : :
1471 [ + - ]: 233 : if (++m_nchk == m_nelem.size()) { // all worker arrays have checkpointed
1472 : 233 : m_nchk = 0;
1473 : : #ifndef HAS_EXAM2M
1474 : : const auto benchmark = g_inputdeck.get< tag::cmd, tag::benchmark >();
1475 : : if (!benchmark) {
1476 : : const auto& restart = g_inputdeck.get< tag::cmd, tag::io, tag::restart >();
1477 : : CkCallback res( CkIndex_Transporter::resume(), thisProxy );
1478 : : CkStartCheckpoint( restart.c_str(), res );
1479 : : } else {
1480 : : resume();
1481 : : }
1482 : : #else
1483 : 233 : resume();
1484 : : #endif
1485 : : }
1486 : 233 : }
1487 : :
1488 : : void
1489 : 233 : Transporter::finish( std::size_t meshid )
1490 : : // *****************************************************************************
1491 : : // Normal finish of time stepping
1492 : : //! \param[in] meshid Mesh id
1493 : : // *****************************************************************************
1494 : : {
1495 : 233 : checkpoint( /* finished = */ 1, meshid );
1496 : 233 : }
1497 : :
1498 : : #include "NoWarning/transporter.def.h"
|
