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