Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/DG.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 DG advances a system of PDEs with the discontinuous Galerkin scheme
9 : : \details DG advances a system of partial differential equations (PDEs) using
10 : : discontinuous Galerkin (DG) finite element (FE) spatial discretization (on
11 : : tetrahedron elements) combined with Runge-Kutta (RK) time stepping.
12 : : \see The documentation in DG.h.
13 : : */
14 : : // *****************************************************************************
15 : :
16 : : #include <algorithm>
17 : : #include <numeric>
18 : : #include <sstream>
19 : :
20 : : #include "DG.hpp"
21 : : #include "Discretization.hpp"
22 : : #include "DGPDE.hpp"
23 : : #include "DiagReducer.hpp"
24 : : #include "DerivedData.hpp"
25 : : #include "ElemDiagnostics.hpp"
26 : : #include "Inciter/InputDeck/InputDeck.hpp"
27 : : #include "Refiner.hpp"
28 : : #include "Limiter.hpp"
29 : : #include "Reorder.hpp"
30 : : #include "Vector.hpp"
31 : : #include "Around.hpp"
32 : : #include "Integrate/Basis.hpp"
33 : : #include "FieldOutput.hpp"
34 : : #include "ChareStateCollector.hpp"
35 : :
36 : : namespace inciter {
37 : :
38 : : extern ctr::InputDeck g_inputdeck;
39 : : extern std::vector< DGPDE > g_dgpde;
40 : :
41 : : //! Runge-Kutta coefficients
42 : : static const std::array< std::array< tk::real, 3 >, 2 >
43 : : rkcoef{{ {{ 0.0, 3.0/4.0, 1.0/3.0 }}, {{ 1.0, 1.0/4.0, 2.0/3.0 }} }};
44 : :
45 : : } // inciter::
46 : :
47 : : extern tk::CProxy_ChareStateCollector stateProxy;
48 : :
49 : : using inciter::DG;
50 : :
51 : 853 : DG::DG( const CProxy_Discretization& disc,
52 : : const CProxy_Ghosts& ghostsproxy,
53 : : const std::map< int, std::vector< std::size_t > >& bface,
54 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
55 : 853 : const std::vector< std::size_t >& triinpoel ) :
56 : : m_disc( disc ),
57 : : m_ghosts( ghostsproxy ),
58 : : m_ndof_NodalExtrm( 3 ), // for the first order derivatives in 3 directions
59 : : m_nsol( 0 ),
60 : : m_ninitsol( 0 ),
61 : : m_nlim( 0 ),
62 : : m_nnod( 0 ),
63 : : m_nrefine( 0 ),
64 : : m_nsmooth( 0 ),
65 : : m_nreco( 0 ),
66 : : m_nnodalExtrema( 0 ),
67 [ + - ]: 853 : m_u( Disc()->Inpoel().size()/4,
68 : 853 : g_inputdeck.get< tag::rdof >()*
69 : 853 : g_inputdeck.get< tag::ncomp >() ),
70 : : m_un( m_u.nunk(), m_u.nprop() ),
71 : 853 : m_p( m_u.nunk(), g_inputdeck.get< tag::rdof >()*
72 [ + - ][ + - ]: 853 : g_dgpde[Disc()->MeshId()].nprim() ),
73 : : m_lhs( m_u.nunk(),
74 : 853 : g_inputdeck.get< tag::ndof >()*
75 : 853 : g_inputdeck.get< tag::ncomp >() ),
76 : : m_rhs( m_u.nunk(), m_lhs.nprop() ),
77 : : m_mtInv(
78 : : tk::invMassMatTaylorRefEl(g_inputdeck.get< tag::rdof >()) ),
79 : : m_uNodalExtrm(),
80 : : m_pNodalExtrm(),
81 : : m_uNodalExtrmc(),
82 : : m_pNodalExtrmc(),
83 [ + - ]: 853 : m_npoin( Disc()->Coord()[0].size() ),
84 : : m_diag(),
85 : : m_stage( 0 ),
86 : : m_ndof(),
87 : : m_numEqDof(),
88 : : m_uc(),
89 : : m_pc(),
90 : : m_ndofc(),
91 : : m_initial( 1 ),
92 : : m_uElemfields( m_u.nunk(),
93 : : g_inputdeck.get< tag::ncomp >() ),
94 : : m_pElemfields( m_u.nunk(),
95 : 853 : m_p.nprop() / g_inputdeck.get< tag::rdof >() ),
96 : : m_uNodefields( m_npoin,
97 : : g_inputdeck.get< tag::ncomp >() ),
98 : : m_pNodefields( m_npoin,
99 : 853 : m_p.nprop() / g_inputdeck.get< tag::rdof >() ),
100 : : m_uNodefieldsc(),
101 : : m_pNodefieldsc(),
102 : : m_outmesh(),
103 : : m_boxelems(),
104 [ + - ][ + - ]: 4265 : m_shockmarker(m_u.nunk(), 1)
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
105 : : // *****************************************************************************
106 : : // Constructor
107 : : //! \param[in] disc Discretization proxy
108 : : //! \param[in] bface Boundary-faces mapped to side set ids
109 : : //! \param[in] triinpoel Boundary-face connectivity
110 : : // *****************************************************************************
111 : : {
112 [ + + ]: 853 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
113 [ + + ]: 815 : g_inputdeck.get< tag::cmd, tag::quiescence >())
114 [ + - ][ + - ]: 1016 : stateProxy.ckLocalBranch()->insert( "DG", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
[ - + ]
115 : : "DG" );
116 : :
117 : : // assign number of dofs for each equation in all pde systems
118 [ + - ][ + - ]: 853 : g_dgpde[Disc()->MeshId()].numEquationDofs(m_numEqDof);
119 : :
120 : : // Allocate storage for the vector of nodal extrema
121 [ + - ][ + - ]: 853 : m_uNodalExtrm.resize( Disc()->Bid().size(),
122 : 0 : std::vector<tk::real>( 2 * m_ndof_NodalExtrm *
123 [ + - ]: 853 : g_inputdeck.get< tag::ncomp >() ) );
124 [ + - ][ + - ]: 853 : m_pNodalExtrm.resize( Disc()->Bid().size(),
125 : 0 : std::vector<tk::real>( 2 * m_ndof_NodalExtrm *
126 [ + - ]: 853 : m_p.nprop() / g_inputdeck.get< tag::rdof >() ) );
127 : :
128 : : // Initialization for the buffer vector of nodal extrema
129 [ + - ]: 853 : resizeNodalExtremac();
130 : :
131 : 853 : usesAtSync = true; // enable migration at AtSync
132 : :
133 : : // Enable SDAG wait for initially building the solution vector and limiting
134 [ + - ]: 853 : if (m_initial) {
135 [ + - ][ + - ]: 853 : thisProxy[ thisIndex ].wait4sol();
136 [ + - ][ + - ]: 853 : thisProxy[ thisIndex ].wait4refine();
137 [ + - ][ + - ]: 853 : thisProxy[ thisIndex ].wait4smooth();
138 [ + - ][ + - ]: 853 : thisProxy[ thisIndex ].wait4lim();
139 [ + - ][ + - ]: 853 : thisProxy[ thisIndex ].wait4nod();
140 [ + - ][ + - ]: 853 : thisProxy[ thisIndex ].wait4reco();
141 [ + - ][ + - ]: 1706 : thisProxy[ thisIndex ].wait4nodalExtrema();
142 : : }
143 : :
144 [ + - ][ + - ]: 1706 : m_ghosts[thisIndex].insert(m_disc, bface, triinpoel, m_u.nunk(),
145 [ + - ][ + - ]: 2559 : CkCallback(CkIndex_DG::resizeSolVectors(), thisProxy[thisIndex]));
[ - + ][ - + ]
[ - - ][ - - ]
146 : :
147 : : // global-sync to call doneInserting on m_ghosts
148 [ + - ]: 853 : auto meshid = Disc()->MeshId();
149 [ + - ]: 853 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
150 [ + - ][ - - ]: 853 : CkCallback(CkReductionTarget(Transporter,doneInsertingGhosts),
151 [ + - ]: 853 : Disc()->Tr()) );
152 : 853 : }
153 : :
154 : : void
155 : 587 : DG::registerReducers()
156 : : // *****************************************************************************
157 : : // Configure Charm++ reduction types
158 : : //! \details Since this is a [initnode] routine, the runtime system executes the
159 : : //! routine exactly once on every logical node early on in the Charm++ init
160 : : //! sequence. Must be static as it is called without an object. See also:
161 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
162 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
163 : : // *****************************************************************************
164 : : {
165 : 587 : ElemDiagnostics::registerReducers();
166 : 587 : }
167 : :
168 : : void
169 : 20187 : DG::ResumeFromSync()
170 : : // *****************************************************************************
171 : : // Return from migration
172 : : //! \details This is called when load balancing (LB) completes. The presence of
173 : : //! this function does not affect whether or not we block on LB.
174 : : // *****************************************************************************
175 : : {
176 [ - + ][ - - ]: 20187 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
177 : :
178 [ + - ]: 20187 : if (!g_inputdeck.get< tag::cmd, tag::nonblocking >()) next();
179 : 20187 : }
180 : :
181 : : void
182 : 853 : DG::resizeSolVectors()
183 : : // *****************************************************************************
184 : : // Resize solution vectors after extension due to Ghosts and continue with setup
185 : : // *****************************************************************************
186 : : {
187 : : // Resize solution vectors, lhs and rhs by the number of ghost tets
188 : 853 : m_u.resize( myGhosts()->m_nunk );
189 : 853 : m_un.resize( myGhosts()->m_nunk );
190 : 853 : m_p.resize( myGhosts()->m_nunk );
191 : 853 : m_lhs.resize( myGhosts()->m_nunk );
192 : 853 : m_rhs.resize( myGhosts()->m_nunk );
193 : :
194 : : // Size communication buffer for solution and number of degrees of freedom
195 [ + + ]: 3412 : for (auto& n : m_ndofc) n.resize( myGhosts()->m_bid.size() );
196 [ + + ]: 3412 : for (auto& u : m_uc) u.resize( myGhosts()->m_bid.size() );
197 [ + + ]: 3412 : for (auto& p : m_pc) p.resize( myGhosts()->m_bid.size() );
198 : :
199 : : // Initialize number of degrees of freedom in mesh elements
200 : 853 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
201 [ + + ]: 853 : if( pref )
202 : : {
203 : 134 : const auto ndofmax = g_inputdeck.get< tag::pref, tag::ndofmax >();
204 : 134 : m_ndof.resize( myGhosts()->m_nunk, ndofmax );
205 : : }
206 : : else
207 : : {
208 : 719 : const auto ndof = g_inputdeck.get< tag::ndof >();
209 : 719 : m_ndof.resize( myGhosts()->m_nunk, ndof );
210 : : }
211 : :
212 : : // Ensure that we also have all the geometry and connectivity data
213 : : // (including those of ghosts)
214 : : Assert( myGhosts()->m_geoElem.nunk() == m_u.nunk(),
215 : : "GeoElem unknowns size mismatch" );
216 : :
217 : : // Signal the runtime system that all workers have received their adjacency
218 : 853 : std::vector< std::size_t > meshdata{ myGhosts()->m_initial, Disc()->MeshId() };
219 : 853 : contribute( meshdata, CkReduction::sum_ulong,
220 [ + - ][ + - ]: 2559 : CkCallback(CkReductionTarget(Transporter,comfinal), Disc()->Tr()) );
[ + - ][ - - ]
221 : 853 : }
222 : :
223 : : void
224 : 853 : DG::setup()
225 : : // *****************************************************************************
226 : : // Set initial conditions, generate lhs, output mesh
227 : : // *****************************************************************************
228 : : {
229 [ + + ]: 853 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
230 [ + + ]: 815 : g_inputdeck.get< tag::cmd, tag::quiescence >())
231 [ + - ][ + - ]: 1016 : stateProxy.ckLocalBranch()->insert( "DG", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
[ - + ]
232 : : "setup" );
233 : :
234 : 853 : auto d = Disc();
235 : :
236 : : // Basic error checking on sizes of element geometry data and connectivity
237 : : Assert( myGhosts()->m_geoElem.nunk() == m_lhs.nunk(),
238 : : "Size mismatch in DG::setup()" );
239 : :
240 : : // Compute left-hand side of discrete PDEs
241 : 853 : lhs();
242 : :
243 : : // Determine elements inside user-defined IC box
244 : 853 : g_dgpde[d->MeshId()].IcBoxElems( myGhosts()->m_geoElem,
245 : 853 : myGhosts()->m_fd.Esuel().size()/4, m_boxelems );
246 : :
247 : : // Compute volume of user-defined box IC
248 [ + - ][ - + ]: 853 : d->boxvol( {}, {}, 0 ); // punt for now
249 : :
250 : : // Query time history field output labels from all PDEs integrated
251 : : const auto& hist_points = g_inputdeck.get< tag::history_output, tag::point >();
252 [ - + ]: 853 : if (!hist_points.empty()) {
253 : 0 : std::vector< std::string > histnames;
254 [ - - ]: 0 : auto n = g_dgpde[d->MeshId()].histNames();
255 [ - - ]: 0 : histnames.insert( end(histnames), begin(n), end(n) );
256 [ - - ]: 0 : d->histheader( std::move(histnames) );
257 : : }
258 : 853 : }
259 : :
260 : : void
261 : 853 : DG::box( tk::real v, const std::vector< tk::real >& )
262 : : // *****************************************************************************
263 : : // Receive total box IC volume and set conditions in box
264 : : //! \param[in] v Total volume within user-specified box
265 : : // *****************************************************************************
266 : : {
267 : 853 : auto d = Disc();
268 : :
269 : : // Store user-defined box IC volume
270 : 853 : d->Boxvol() = v;
271 : :
272 : : // Set initial conditions for all PDEs
273 : 1706 : g_dgpde[d->MeshId()].initialize( m_lhs, myGhosts()->m_inpoel,
274 : 853 : myGhosts()->m_coord, m_boxelems, d->ElemBlockId(), m_u, d->T(),
275 : 853 : myGhosts()->m_fd.Esuel().size()/4 );
276 : 853 : g_dgpde[d->MeshId()].updatePrimitives( m_u, m_lhs, myGhosts()->m_geoElem, m_p,
277 : 853 : myGhosts()->m_fd.Esuel().size()/4 );
278 : :
279 : : m_un = m_u;
280 : :
281 : : // Output initial conditions to file (regardless of whether it was requested)
282 [ + - ][ + - ]: 2559 : startFieldOutput( CkCallback(CkIndex_DG::start(), thisProxy[thisIndex]) );
[ - + ][ - - ]
283 : 853 : }
284 : :
285 : : void
286 : 853 : DG::start()
287 : : // *****************************************************************************
288 : : // Start time stepping
289 : : // *****************************************************************************
290 : : {
291 : : // Free memory storing output mesh
292 : 853 : m_outmesh.destroy();
293 : :
294 : : // Start timer measuring time stepping wall clock time
295 : 853 : Disc()->Timer().zero();
296 : : // Zero grind-timer
297 : 853 : Disc()->grindZero();
298 : : // Start time stepping by computing the size of the next time step)
299 : 853 : next();
300 : 853 : }
301 : :
302 : : void
303 : 25013 : DG::startFieldOutput( CkCallback c )
304 : : // *****************************************************************************
305 : : // Start preparing fields for output to file
306 : : //! \param[in] c Function to continue with after the write
307 : : // *****************************************************************************
308 : : {
309 : : // No field output in benchmark mode or if field output frequency not hit
310 [ + + ][ + + ]: 25013 : if (g_inputdeck.get< tag::cmd, tag::benchmark >() || !fieldOutput()) {
311 : :
312 : 22837 : c.send();
313 : :
314 : : } else {
315 : :
316 : : // Optionally refine mesh for field output
317 : 2176 : auto d = Disc();
318 : :
319 [ + + ]: 2176 : if (refinedOutput()) {
320 : :
321 : 33 : const auto& tr = tk::remap( myGhosts()->m_fd.Triinpoel(), d->Gid() );
322 [ + - ][ + - ]: 66 : d->Ref()->outref( myGhosts()->m_fd.Bface(), {}, tr, c );
[ + - ][ + - ]
[ - - ]
323 : :
324 : : } else {
325 : :
326 : : // cut off ghosts from mesh connectivity and coordinates
327 : 2143 : const auto& tr = tk::remap( myGhosts()->m_fd.Triinpoel(), d->Gid() );
328 [ + - ][ + - ]: 4286 : extractFieldOutput( {}, d->Chunk(), d->Coord(), {}, {},
[ + + ][ - - ]
329 [ + - ]: 2143 : d->NodeCommMap(), myGhosts()->m_fd.Bface(), {}, tr, c );
330 : :
331 : : }
332 : :
333 : : }
334 : 25013 : }
335 : :
336 : : void
337 : 72480 : DG::next()
338 : : // *****************************************************************************
339 : : // Advance equations to next time step
340 : : // *****************************************************************************
341 : : {
342 : 72480 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
343 : :
344 : 72480 : auto d = Disc();
345 : :
346 [ + + ][ + + ]: 72480 : if (pref && m_stage == 0 && d->T() > 0)
[ + + ]
347 : 3272 : g_dgpde[d->MeshId()].eval_ndof( myGhosts()->m_nunk, myGhosts()->m_coord,
348 : 1636 : myGhosts()->m_inpoel,
349 : 1636 : myGhosts()->m_fd, m_u, m_p,
350 : : g_inputdeck.get< tag::pref, tag::indicator >(),
351 : : g_inputdeck.get< tag::ndof >(),
352 : : g_inputdeck.get< tag::pref, tag::ndofmax >(),
353 : : g_inputdeck.get< tag::pref, tag::tolref >(),
354 : 1636 : m_ndof );
355 : :
356 : : // communicate solution ghost data (if any)
357 [ + + ]: 72480 : if (myGhosts()->m_sendGhost.empty())
358 : 4335 : comsol_complete();
359 : : else
360 [ + + ]: 636060 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
361 [ + - ]: 499770 : std::vector< std::size_t > tetid( ghostdata.size() );
362 [ + - ][ + - ]: 999540 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
[ - - ]
363 [ + - ]: 999540 : prim( ghostdata.size() );
364 : : std::vector< std::size_t > ndof;
365 : : std::size_t j = 0;
366 [ + + ]: 10020420 : for(const auto& i : ghostdata) {
367 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4,
368 : : "Sending solution ghost data" );
369 [ + - ]: 9520650 : tetid[j] = i;
370 [ + - ][ - + ]: 9520650 : u[j] = m_u[i];
371 [ + - ][ - + ]: 9520650 : prim[j] = m_p[i];
372 [ + + ][ + + ]: 9520650 : if (pref && m_stage == 0) ndof.push_back( m_ndof[i] );
[ + - ]
373 : 9520650 : ++j;
374 : : }
375 [ + - ][ + - ]: 999540 : thisProxy[ cid ].comsol( thisIndex, m_stage, tetid, u, prim, ndof );
[ + + ][ - - ]
376 : : }
377 : :
378 : 72480 : ownsol_complete();
379 : 72480 : }
380 : :
381 : : void
382 : 499770 : DG::comsol( int fromch,
383 : : std::size_t fromstage,
384 : : const std::vector< std::size_t >& tetid,
385 : : const std::vector< std::vector< tk::real > >& u,
386 : : const std::vector< std::vector< tk::real > >& prim,
387 : : const std::vector< std::size_t >& ndof )
388 : : // *****************************************************************************
389 : : // Receive chare-boundary solution ghost data from neighboring chares
390 : : //! \param[in] fromch Sender chare id
391 : : //! \param[in] fromstage Sender chare time step stage
392 : : //! \param[in] tetid Ghost tet ids we receive solution data for
393 : : //! \param[in] u Solution ghost data
394 : : //! \param[in] prim Primitive variables in ghost cells
395 : : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
396 : : //! \details This function receives contributions to the unlimited solution
397 : : //! from fellow chares.
398 : : // *****************************************************************************
399 : : {
400 : : Assert( u.size() == tetid.size(), "Size mismatch in DG::comsol()" );
401 : : Assert( prim.size() == tetid.size(), "Size mismatch in DG::comsol()" );
402 : :
403 : 499770 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
404 : :
405 : 499770 : if (pref && fromstage == 0)
406 : : Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comsol()" );
407 : :
408 : : // Find local-to-ghost tet id map for sender chare
409 : 499770 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
410 : :
411 [ + + ]: 10020420 : for (std::size_t i=0; i<tetid.size(); ++i) {
412 : 9520650 : auto j = tk::cref_find( n, tetid[i] );
413 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
414 : : "Receiving solution non-ghost data" );
415 : 9520650 : auto b = tk::cref_find( myGhosts()->m_bid, j );
416 : : Assert( b < m_uc[0].size(), "Indexing out of bounds" );
417 : 9520650 : m_uc[0][b] = u[i];
418 : 9520650 : m_pc[0][b] = prim[i];
419 [ + + ]: 9520650 : if (pref && fromstage == 0) {
420 : : Assert( b < m_ndofc[0].size(), "Indexing out of bounds" );
421 : 395110 : m_ndofc[0][b] = ndof[i];
422 : : }
423 : : }
424 : :
425 : : // if we have received all solution ghost contributions from neighboring
426 : : // chares (chares we communicate along chare-boundary faces with), and
427 : : // contributed our solution to these neighbors, proceed to reconstructions
428 [ + + ]: 499770 : if (++m_nsol == myGhosts()->m_sendGhost.size()) {
429 : 68145 : m_nsol = 0;
430 : 68145 : comsol_complete();
431 : : }
432 : 499770 : }
433 : :
434 : : void
435 : 2176 : DG::extractFieldOutput(
436 : : const std::vector< std::size_t >& /*ginpoel*/,
437 : : const tk::UnsMesh::Chunk& chunk,
438 : : const tk::UnsMesh::Coords& coord,
439 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /*addedNodes*/,
440 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
441 : : const tk::NodeCommMap& nodeCommMap,
442 : : const std::map< int, std::vector< std::size_t > >& bface,
443 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
444 : : const std::vector< std::size_t >& triinpoel,
445 : : CkCallback c )
446 : : // *****************************************************************************
447 : : // Extract field output going to file
448 : : //! \param[in] chunk Field-output mesh chunk (connectivity and global<->local
449 : : //! id maps)
450 : : //! \param[in] coord Field-output mesh node coordinates
451 : : //! \param[in] addedTets Field-output mesh cells and their parents (local ids)
452 : : //! \param[in] nodeCommMap Field-output mesh node communication map
453 : : //! \param[in] bface Field-output meshndary-faces mapped to side set ids
454 : : //! \param[in] triinpoel Field-output mesh boundary-face connectivity
455 : : //! \param[in] c Function to continue with after the write
456 : : // *****************************************************************************
457 : : {
458 : : m_outmesh.chunk = chunk;
459 : : m_outmesh.coord = coord;
460 : 2176 : m_outmesh.triinpoel = triinpoel;
461 : : m_outmesh.bface = bface;
462 : : m_outmesh.nodeCommMap = nodeCommMap;
463 : :
464 : : const auto& inpoel = std::get< 0 >( chunk );
465 : :
466 : : // Evaluate element solution on incoming mesh
467 : 2176 : evalSolution( *Disc(), inpoel, coord, addedTets, m_ndof, m_u, m_p,
468 : 2176 : m_uElemfields, m_pElemfields, m_uNodefields, m_pNodefields );
469 : :
470 : : // Send node fields contributions to neighbor chares
471 [ + + ]: 2176 : if (nodeCommMap.empty())
472 : 180 : comnodeout_complete();
473 : : else {
474 : : const auto& lid = std::get< 2 >( chunk );
475 : 3992 : auto esup = tk::genEsup( inpoel, 4 );
476 [ + + ]: 15326 : for(const auto& [ch,nodes] : nodeCommMap) {
477 : : // Pack node field data in chare boundary nodes
478 : : std::vector< std::vector< tk::real > >
479 [ + - ][ + - ]: 39990 : lu( m_uNodefields.nprop(), std::vector< tk::real >( nodes.size() ) );
[ + - ]
480 : : std::vector< std::vector< tk::real > >
481 [ + - ][ + - ]: 39990 : lp( m_pNodefields.nprop(), std::vector< tk::real >( nodes.size() ) );
482 [ + + ]: 62724 : for (std::size_t f=0; f<m_uNodefields.nprop(); ++f) {
483 : : std::size_t j = 0;
484 [ + + ]: 509126 : for (auto g : nodes)
485 : 459732 : lu[f][j++] = m_uNodefields(tk::cref_find(lid,g),f);
486 : : }
487 [ + + ]: 21242 : for (std::size_t f=0; f<m_pNodefields.nprop(); ++f) {
488 : : std::size_t j = 0;
489 [ + + ]: 88936 : for (auto g : nodes)
490 : 81024 : lp[f][j++] = m_pNodefields(tk::cref_find(lid,g),f);
491 : : }
492 : : // Pack (partial) number of elements surrounding chare boundary nodes
493 [ + - ]: 13330 : std::vector< std::size_t > nesup( nodes.size() );
494 : : std::size_t j = 0;
495 [ + + ]: 141372 : for (auto g : nodes) {
496 : 128042 : auto i = tk::cref_find( lid, g );
497 : 128042 : nesup[j++] = esup.second[i+1] - esup.second[i];
498 : : }
499 [ + - ][ + - ]: 39990 : thisProxy[ch].comnodeout(
[ + - ][ - - ]
500 [ + - ][ - + ]: 26660 : std::vector<std::size_t>(begin(nodes),end(nodes)), nesup, lu, lp );
501 : : }
502 : : }
503 : :
504 [ + - ]: 2176 : ownnod_complete( c, addedTets );
505 : 2176 : }
506 : :
507 : : void
508 : 853 : DG::lhs()
509 : : // *****************************************************************************
510 : : // Compute left-hand side of discrete transport equations
511 : : // *****************************************************************************
512 : : {
513 : 853 : g_dgpde[Disc()->MeshId()].lhs( myGhosts()->m_geoElem, m_lhs );
514 : :
515 [ - + ]: 853 : if (!m_initial) stage();
516 : 853 : }
517 : :
518 : 72480 : void DG::refine()
519 : : // *****************************************************************************
520 : : // Add the protective layer for ndof refinement
521 : : // *****************************************************************************
522 : : {
523 : 72480 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
524 : :
525 : : // Combine own and communicated contributions of unreconstructed solution and
526 : : // degrees of freedom in cells (if p-adaptive)
527 [ + + ]: 9665610 : for (const auto& b : myGhosts()->m_bid) {
528 : : Assert( m_uc[0][b.second].size() == m_u.nprop(), "ncomp size mismatch" );
529 : : Assert( m_pc[0][b.second].size() == m_p.nprop(), "ncomp size mismatch" );
530 [ + + ]: 149542125 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
531 : 140021475 : m_u(b.first,c) = m_uc[0][b.second][c];
532 : : }
533 [ + + ]: 30332715 : for (std::size_t c=0; c<m_p.nprop(); ++c) {
534 : 20812065 : m_p(b.first,c) = m_pc[0][b.second][c];
535 : : }
536 [ + + ][ + + ]: 9520650 : if (pref && m_stage == 0) {
537 : 395110 : m_ndof[ b.first ] = m_ndofc[0][ b.second ];
538 : : }
539 : : }
540 : :
541 [ + + ][ + + ]: 72480 : if (pref && m_stage==0) refine_ndof();
542 : :
543 [ + + ]: 72480 : if (myGhosts()->m_sendGhost.empty())
544 : 4335 : comrefine_complete();
545 : : else
546 [ + + ]: 636060 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
547 [ + - ]: 499770 : std::vector< std::size_t > tetid( ghostdata.size() );
548 [ + - ][ + - ]: 999540 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
[ - - ]
549 [ + - ][ + - ]: 999540 : prim( ghostdata.size() );
550 [ + - ]: 499770 : std::vector< std::size_t > ndof( ghostdata.size() );
551 : : std::size_t j = 0;
552 [ + + ]: 10020420 : for(const auto& i : ghostdata) {
553 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4, "Sending refined ndof "
554 : : "data" );
555 [ + + ]: 9520650 : tetid[j] = i;
556 [ + + ][ + + ]: 9520650 : if (pref && m_stage == 0) ndof[j] = m_ndof[i];
557 : 9520650 : ++j;
558 : : }
559 [ + - ][ + - ]: 999540 : thisProxy[ cid ].comrefine( thisIndex, tetid, ndof );
[ + - ][ - - ]
560 : : }
561 : :
562 : 72480 : ownrefine_complete();
563 : 72480 : }
564 : :
565 : : void
566 : 499770 : DG::comrefine( int fromch,
567 : : const std::vector< std::size_t >& tetid,
568 : : const std::vector< std::size_t >& ndof )
569 : : // *****************************************************************************
570 : : // Receive chare-boundary ghost data from neighboring chares
571 : : //! \param[in] fromch Sender chare id
572 : : //! \param[in] tetid Ghost tet ids we receive solution data for
573 : : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
574 : : //! \details This function receives contributions to the refined ndof data
575 : : //! from fellow chares.
576 : : // *****************************************************************************
577 : : {
578 : 499770 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
579 : :
580 : : if (pref && m_stage == 0)
581 : : Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comrefine()" );
582 : :
583 : : // Find local-to-ghost tet id map for sender chare
584 : 499770 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
585 : :
586 [ + + ]: 10020420 : for (std::size_t i=0; i<tetid.size(); ++i) {
587 : 9520650 : auto j = tk::cref_find( n, tetid[i] );
588 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
589 : : "Receiving solution non-ghost data" );
590 : 9520650 : auto b = tk::cref_find( myGhosts()->m_bid, j );
591 [ + + ][ + + ]: 9520650 : if (pref && m_stage == 0) {
592 : : Assert( b < m_ndofc[1].size(), "Indexing out of bounds" );
593 : 395110 : m_ndofc[1][b] = ndof[i];
594 : : }
595 : : }
596 : :
597 : : // if we have received all solution ghost contributions from neighboring
598 : : // chares (chares we communicate along chare-boundary faces with), and
599 : : // contributed our solution to these neighbors, proceed to limiting
600 [ + + ]: 499770 : if (++m_nrefine == myGhosts()->m_sendGhost.size()) {
601 : 68145 : m_nrefine = 0;
602 : 68145 : comrefine_complete();
603 : : }
604 : 499770 : }
605 : :
606 : : void
607 : 72480 : DG::smooth()
608 : : // *****************************************************************************
609 : : // Smooth the refined ndof distribution
610 : : // *****************************************************************************
611 : : {
612 : 72480 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
613 : :
614 [ + + ]: 9665610 : for (const auto& b : myGhosts()->m_bid) {
615 [ + + ][ + + ]: 9520650 : if (pref && m_stage == 0)
616 : 395110 : m_ndof[ b.first ] = m_ndofc[1][ b.second ];
617 : : }
618 : :
619 [ + + ][ + + ]: 72480 : if (pref && m_stage==0) smooth_ndof();
620 : :
621 [ + + ]: 72480 : if (myGhosts()->m_sendGhost.empty())
622 : 4335 : comsmooth_complete();
623 : : else
624 [ + + ]: 636060 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
625 : 499770 : std::vector< std::size_t > tetid( ghostdata.size() );
626 : : std::vector< std::size_t > ndof;
627 : : std::size_t j = 0;
628 [ + + ]: 10020420 : for(const auto& i : ghostdata) {
629 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4, "Sending ndof data" );
630 [ + + ]: 9520650 : tetid[j] = i;
631 [ + + ][ + + ]: 9520650 : if (pref && m_stage == 0) ndof.push_back( m_ndof[i] );
[ + - ]
632 : 9520650 : ++j;
633 : : }
634 [ + - ][ + - ]: 999540 : thisProxy[ cid ].comsmooth( thisIndex, tetid, ndof );
[ + + ][ - - ]
635 : : }
636 : :
637 : 72480 : ownsmooth_complete();
638 : 72480 : }
639 : :
640 : : void
641 : 499770 : DG::comsmooth( int fromch,
642 : : const std::vector< std::size_t >& tetid,
643 : : const std::vector< std::size_t >& ndof )
644 : : // *****************************************************************************
645 : : // Receive chare-boundary ghost data from neighboring chares
646 : : //! \param[in] fromch Sender chare id
647 : : //! \param[in] tetid Ghost tet ids we receive solution data for
648 : : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
649 : : //! \details This function receives contributions to the smoothed ndof data
650 : : //! from fellow chares.
651 : : // *****************************************************************************
652 : : {
653 : 499770 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
654 : :
655 : : if (pref && m_stage == 0)
656 : : Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comsmooth()" );
657 : :
658 : 499770 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
659 : :
660 [ + + ]: 10020420 : for (std::size_t i=0; i<tetid.size(); ++i) {
661 : 9520650 : auto j = tk::cref_find( n, tetid[i] );
662 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4, "Receiving ndof data" );
663 : 9520650 : auto b = tk::cref_find( myGhosts()->m_bid, j );
664 [ + + ][ + + ]: 9520650 : if (pref && m_stage == 0) {
665 : : Assert( b < m_ndofc[2].size(), "Indexing out of bounds" );
666 : 395110 : m_ndofc[2][b] = ndof[i];
667 : : }
668 : : }
669 : :
670 [ + + ]: 499770 : if (++m_nsmooth == myGhosts()->m_sendGhost.size()) {
671 : 68145 : m_nsmooth = 0;
672 : 68145 : comsmooth_complete();
673 : : }
674 : 499770 : }
675 : :
676 : : void
677 : 72480 : DG::reco()
678 : : // *****************************************************************************
679 : : // Compute reconstructions
680 : : // *****************************************************************************
681 : : {
682 : 72480 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
683 : 72480 : const auto rdof = g_inputdeck.get< tag::rdof >();
684 : :
685 : : // Combine own and communicated contributions of unreconstructed solution and
686 : : // degrees of freedom in cells (if p-adaptive)
687 [ + + ]: 9665610 : for (const auto& b : myGhosts()->m_bid) {
688 [ + + ][ + + ]: 9520650 : if (pref && m_stage == 0) {
689 : 395110 : m_ndof[ b.first ] = m_ndofc[2][ b.second ];
690 : : }
691 : : }
692 : :
693 : 72480 : auto d = Disc();
694 [ + + ]: 72480 : if (rdof > 1)
695 : : // Reconstruct second-order solution and primitive quantities
696 : 89280 : g_dgpde[d->MeshId()].reconstruct( d->T(), myGhosts()->m_geoFace,
697 : 44640 : myGhosts()->m_geoElem,
698 : 44640 : myGhosts()->m_fd, myGhosts()->m_esup, myGhosts()->m_inpoel,
699 : 44640 : myGhosts()->m_coord, m_u, m_p );
700 : :
701 : : // Send reconstructed solution to neighboring chares
702 [ + + ]: 72480 : if (myGhosts()->m_sendGhost.empty())
703 : 4335 : comreco_complete();
704 : : else
705 [ + + ]: 636060 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
706 [ + - ]: 499770 : std::vector< std::size_t > tetid( ghostdata.size() );
707 [ + - ][ + - ]: 999540 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
[ - - ]
708 [ + - ]: 499770 : prim( ghostdata.size() );
709 : : std::size_t j = 0;
710 [ + + ]: 10020420 : for(const auto& i : ghostdata) {
711 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4, "Sending reconstructed ghost "
712 : : "data" );
713 [ + - ]: 9520650 : tetid[j] = i;
714 [ + - ][ - + ]: 9520650 : u[j] = m_u[i];
715 [ + - ][ - + ]: 9520650 : prim[j] = m_p[i];
716 : 9520650 : ++j;
717 : : }
718 [ + - ][ + - ]: 999540 : thisProxy[ cid ].comreco( thisIndex, tetid, u, prim );
719 : : }
720 : :
721 : 72480 : ownreco_complete();
722 : 72480 : }
723 : :
724 : : void
725 : 499770 : DG::comreco( int fromch,
726 : : const std::vector< std::size_t >& tetid,
727 : : const std::vector< std::vector< tk::real > >& u,
728 : : const std::vector< std::vector< tk::real > >& prim )
729 : : // *****************************************************************************
730 : : // Receive chare-boundary reconstructed ghost data from neighboring chares
731 : : //! \param[in] fromch Sender chare id
732 : : //! \param[in] tetid Ghost tet ids we receive solution data for
733 : : //! \param[in] u Reconstructed high-order solution
734 : : //! \param[in] prim Limited high-order primitive quantities
735 : : //! \details This function receives contributions to the reconstructed solution
736 : : //! from fellow chares.
737 : : // *****************************************************************************
738 : : {
739 : : Assert( u.size() == tetid.size(), "Size mismatch in DG::comreco()" );
740 : : Assert( prim.size() == tetid.size(), "Size mismatch in DG::comreco()" );
741 : :
742 : : // Find local-to-ghost tet id map for sender chare
743 : 499770 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
744 : :
745 [ + + ]: 10020420 : for (std::size_t i=0; i<tetid.size(); ++i) {
746 : 9520650 : auto j = tk::cref_find( n, tetid[i] );
747 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
748 : : "Receiving solution non-ghost data" );
749 : 9520650 : auto b = tk::cref_find( myGhosts()->m_bid, j );
750 : : Assert( b < m_uc[1].size(), "Indexing out of bounds" );
751 : : Assert( b < m_pc[1].size(), "Indexing out of bounds" );
752 : 9520650 : m_uc[1][b] = u[i];
753 : 9520650 : m_pc[1][b] = prim[i];
754 : : }
755 : :
756 : : // if we have received all solution ghost contributions from neighboring
757 : : // chares (chares we communicate along chare-boundary faces with), and
758 : : // contributed our solution to these neighbors, proceed to limiting
759 [ + + ]: 499770 : if (++m_nreco == myGhosts()->m_sendGhost.size()) {
760 : 68145 : m_nreco = 0;
761 : 68145 : comreco_complete();
762 : : }
763 : 499770 : }
764 : :
765 : : void
766 : 72480 : DG::nodalExtrema()
767 : : // *****************************************************************************
768 : : // Compute nodal extrema at chare-boundary nodes. Extrema at internal nodes
769 : : // are calculated in limiter function.
770 : : // *****************************************************************************
771 : : {
772 : 72480 : auto d = Disc();
773 : 72480 : auto gid = d->Gid();
774 : : auto bid = d->Bid();
775 : 72480 : const auto rdof = g_inputdeck.get< tag::rdof >();
776 : 72480 : const auto ncomp = m_u.nprop() / rdof;
777 : 72480 : const auto nprim = m_p.nprop() / rdof;
778 : :
779 : : // Combine own and communicated contributions of unlimited solution, and
780 : : // if a p-adaptive algorithm is used, degrees of freedom in cells
781 [ + - ][ + + ]: 9665610 : for (const auto& [boundary, localtet] : myGhosts()->m_bid) {
782 : : Assert( m_uc[1][localtet].size() == m_u.nprop(), "ncomp size mismatch" );
783 : : Assert( m_pc[1][localtet].size() == m_p.nprop(), "ncomp size mismatch" );
784 [ + + ]: 149542125 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
785 : 140021475 : m_u(boundary,c) = m_uc[1][localtet][c];
786 : : }
787 [ + + ]: 30332715 : for (std::size_t c=0; c<m_p.nprop(); ++c) {
788 : 20812065 : m_p(boundary,c) = m_pc[1][localtet][c];
789 : : }
790 : : }
791 : :
792 : : // Initialize nodal extrema vector
793 : : auto large = std::numeric_limits< tk::real >::max();
794 [ + + ]: 2360640 : for(std::size_t i = 0; i<bid.size(); i++)
795 : : {
796 [ + + ]: 8291580 : for (std::size_t c=0; c<ncomp; ++c)
797 : : {
798 [ + + ]: 24013680 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
799 : : {
800 : 18010260 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
801 : 18010260 : auto min_mark = max_mark + 1;
802 : 18010260 : m_uNodalExtrm[i][max_mark] = -large;
803 : 18010260 : m_uNodalExtrm[i][min_mark] = large;
804 : : }
805 : : }
806 [ + + ]: 3494370 : for (std::size_t c=0; c<nprim; ++c)
807 : : {
808 [ + + ]: 4824840 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
809 : : {
810 : 3618630 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
811 : 3618630 : auto min_mark = max_mark + 1;
812 : 3618630 : m_pNodalExtrm[i][max_mark] = -large;
813 : 3618630 : m_pNodalExtrm[i][min_mark] = large;
814 : : }
815 : : }
816 : : }
817 : :
818 : : // Evaluate the max/min value for the chare-boundary nodes
819 [ + + ]: 72480 : if(rdof > 4) {
820 [ + - ]: 14970 : evalNodalExtrmRefEl(ncomp, nprim, m_ndof_NodalExtrm, d->bndel(),
821 [ + - ][ + - ]: 7485 : myGhosts()->m_inpoel, gid, bid, m_u, m_p, m_uNodalExtrm, m_pNodalExtrm);
822 : : }
823 : :
824 : : // Communicate extrema at nodes to other chares on chare-boundary
825 [ + + ]: 72480 : if (d->NodeCommMap().empty()) // in serial we are done
826 [ + - ]: 4335 : comnodalExtrema_complete();
827 : : else // send nodal extrema to chare-boundary nodes to fellow chares
828 : : {
829 [ + + ]: 567915 : for (const auto& [c,n] : d->NodeCommMap()) {
830 [ + - ][ + - ]: 999540 : std::vector< std::vector< tk::real > > g1( n.size() ), g2( n.size() );
831 : : std::size_t j = 0;
832 [ + + ]: 4202640 : for (auto i : n)
833 : : {
834 : 3702870 : auto p = tk::cref_find(d->Bid(),i);
835 [ + - ]: 3702870 : g1[ j ] = m_uNodalExtrm[ p ];
836 [ + - ]: 3702870 : g2[ j++ ] = m_pNodalExtrm[ p ];
837 : : }
838 [ + - ][ + - ]: 999540 : thisProxy[c].comnodalExtrema( std::vector<std::size_t>(begin(n),end(n)),
[ + - ][ - + ]
839 : : g1, g2 );
840 : : }
841 : : }
842 [ + - ]: 72480 : ownnodalExtrema_complete();
843 : 72480 : }
844 : :
845 : : void
846 : 499770 : DG::comnodalExtrema( const std::vector< std::size_t >& gid,
847 : : const std::vector< std::vector< tk::real > >& G1,
848 : : const std::vector< std::vector< tk::real > >& G2 )
849 : : // *****************************************************************************
850 : : // Receive contributions to nodal extrema on chare-boundaries
851 : : //! \param[in] gid Global mesh node IDs at which we receive grad contributions
852 : : //! \param[in] G1 Partial contributions of extrema for conservative variables to
853 : : //! chare-boundary nodes
854 : : //! \param[in] G2 Partial contributions of extrema for primitive variables to
855 : : //! chare-boundary nodes
856 : : //! \details This function receives contributions to m_uNodalExtrm/m_pNodalExtrm
857 : : //! , which stores nodal extrems at mesh chare-boundary nodes. While
858 : : //! m_uNodalExtrm/m_pNodalExtrm stores own contributions, m_uNodalExtrmc
859 : : //! /m_pNodalExtrmc collects the neighbor chare contributions during
860 : : //! communication.
861 : : // *****************************************************************************
862 : : {
863 : : Assert( G1.size() == gid.size(), "Size mismatch" );
864 : : Assert( G2.size() == gid.size(), "Size mismatch" );
865 : :
866 : 499770 : const auto rdof = g_inputdeck.get< tag::rdof >();
867 : 499770 : const auto ncomp = m_u.nprop() / rdof;
868 : 499770 : const auto nprim = m_p.nprop() / rdof;
869 : :
870 [ + + ]: 4202640 : for (std::size_t i=0; i<gid.size(); ++i)
871 : : {
872 : : auto& u = m_uNodalExtrmc[gid[i]];
873 : 3702870 : auto& p = m_pNodalExtrmc[gid[i]];
874 [ + + ]: 15381360 : for (std::size_t c=0; c<ncomp; ++c)
875 : : {
876 [ + + ]: 46713960 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
877 : : {
878 : 35035470 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
879 : 35035470 : auto min_mark = max_mark + 1;
880 [ + + ][ + + ]: 35611549 : u[max_mark] = std::max( G1[i][max_mark], u[max_mark] );
881 : 35035470 : u[min_mark] = std::min( G1[i][min_mark], u[min_mark] );
882 : : }
883 : : }
884 [ + + ]: 6489900 : for (std::size_t c=0; c<nprim; ++c)
885 : : {
886 [ + + ]: 11148120 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
887 : : {
888 : 8361090 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
889 : 8361090 : auto min_mark = max_mark + 1;
890 [ - + ][ - + ]: 8361090 : p[max_mark] = std::max( G2[i][max_mark], p[max_mark] );
891 : 8361090 : p[min_mark] = std::min( G2[i][min_mark], p[min_mark] );
892 : : }
893 : : }
894 : : }
895 : :
896 [ + + ]: 499770 : if (++m_nnodalExtrema == Disc()->NodeCommMap().size())
897 : : {
898 : 68145 : m_nnodalExtrema = 0;
899 : 68145 : comnodalExtrema_complete();
900 : : }
901 : 499770 : }
902 : :
903 : 73333 : void DG::resizeNodalExtremac()
904 : : // *****************************************************************************
905 : : // Resize the buffer vector of nodal extrema
906 : : // *****************************************************************************
907 : : {
908 : 73333 : const auto rdof = g_inputdeck.get< tag::rdof >();
909 : 73333 : const auto ncomp = m_u.nprop() / rdof;
910 : 73333 : const auto nprim = m_p.nprop() / rdof;
911 : :
912 : 73333 : auto large = std::numeric_limits< tk::real >::max();
913 [ + - ][ + + ]: 653708 : for (const auto& [c,n] : Disc()->NodeCommMap())
914 : : {
915 [ + + ][ + - ]: 4257712 : for (auto i : n) {
916 : : auto& u = m_uNodalExtrmc[i];
917 : : auto& p = m_pNodalExtrmc[i];
918 [ + - ]: 3750670 : u.resize( 2*m_ndof_NodalExtrm*ncomp, large );
919 [ + - ]: 3750670 : p.resize( 2*m_ndof_NodalExtrm*nprim, large );
920 : :
921 : : // Initialize the minimum nodal extrema
922 [ + + ]: 15002680 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
923 : : {
924 [ + + ]: 46768698 : for(std::size_t k = 0; k < ncomp; k++)
925 : 35516688 : u[2*k*m_ndof_NodalExtrm+2*idof] = -large;
926 [ + + ]: 19725294 : for(std::size_t k = 0; k < nprim; k++)
927 : 8473284 : p[2*k*m_ndof_NodalExtrm+2*idof] = -large;
928 : : }
929 : : }
930 : : }
931 : 73333 : }
932 : :
933 : 7485 : void DG::evalNodalExtrmRefEl(
934 : : const std::size_t ncomp,
935 : : const std::size_t nprim,
936 : : const std::size_t ndof_NodalExtrm,
937 : : const std::vector< std::size_t >& bndel,
938 : : const std::vector< std::size_t >& inpoel,
939 : : const std::vector< std::size_t >& gid,
940 : : const std::unordered_map< std::size_t, std::size_t >& bid,
941 : : const tk::Fields& U,
942 : : const tk::Fields& P,
943 : : std::vector< std::vector<tk::real> >& uNodalExtrm,
944 : : std::vector< std::vector<tk::real> >& pNodalExtrm )
945 : : // *****************************************************************************
946 : : // Compute the nodal extrema of ref el derivatives for chare-boundary nodes
947 : : //! \param[in] ncomp Number of conservative variables
948 : : //! \param[in] nprim Number of primitive variables
949 : : //! \param[in] ndof_NodalExtrm Degree of freedom for nodal extrema
950 : : //! \param[in] bndel List of elements contributing to chare-boundary nodes
951 : : //! \param[in] inpoel Element-node connectivity for element e
952 : : //! \param[in] gid Local->global node id map
953 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
954 : : //! global node ids (key)
955 : : //! \param[in] U Vector of conservative variables
956 : : //! \param[in] P Vector of primitive variables
957 : : //! \param[in,out] uNodalExtrm Chare-boundary nodal extrema for conservative
958 : : //! variables
959 : : //! \param[in,out] pNodalExtrm Chare-boundary nodal extrema for primitive
960 : : //! variables
961 : : // *****************************************************************************
962 : : {
963 : 7485 : const auto rdof = g_inputdeck.get< tag::rdof >();
964 : :
965 [ + + ]: 702165 : for (auto e : bndel)
966 : : {
967 : : // access node IDs
968 : : const std::vector<std::size_t> N
969 [ + - ]: 694680 : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
970 : :
971 : : // Loop over nodes of element e
972 [ + + ]: 3473400 : for(std::size_t ip=0; ip<4; ++ip)
973 : : {
974 : 2778720 : auto i = bid.find( gid[N[ip]] );
975 [ + + ]: 2778720 : if (i != end(bid)) // If ip is the chare boundary point
976 : : {
977 : : // If DG(P2) is applied, find the nodal extrema of the gradients of
978 : : // conservative/primitive variables in the reference element
979 : :
980 : : // Vector used to store the first order derivatives for both
981 : : // conservative and primitive variables
982 [ + - ][ - - ]: 1861455 : std::vector< std::array< tk::real, 3 > > gradc(ncomp, {0.0, 0.0, 0.0});
983 [ + - ][ - - ]: 1861455 : std::vector< std::array< tk::real, 3 > > gradp(ncomp, {0.0, 0.0, 0.0});
984 : :
985 : : // Derivatives of the Dubiner basis
986 : 1861455 : std::array< tk::real, 3 > center {{0.25, 0.25, 0.25}};
987 [ + - ]: 1861455 : auto dBdxi = tk::eval_dBdxi(rdof, center);
988 : :
989 : : // Evaluate the first order derivative
990 [ + + ]: 10741830 : for(std::size_t icomp = 0; icomp < ncomp; icomp++)
991 : : {
992 : 8880375 : auto mark = icomp * rdof;
993 [ + + ]: 35521500 : for(std::size_t idir = 0; idir < 3; idir++)
994 : : {
995 : 26641125 : gradc[icomp][idir] = 0;
996 [ + + ]: 266411250 : for(std::size_t idof = 1; idof < rdof; idof++)
997 : 239770125 : gradc[icomp][idir] += U(e, mark+idof) * dBdxi[idir][idof];
998 : : }
999 : : }
1000 [ - + ]: 1861455 : for(std::size_t icomp = 0; icomp < nprim; icomp++)
1001 : : {
1002 : 0 : auto mark = icomp * rdof;
1003 [ - - ]: 0 : for(std::size_t idir = 0; idir < 3; idir++)
1004 : : {
1005 : 0 : gradp[icomp][idir] = 0;
1006 [ - - ]: 0 : for(std::size_t idof = 1; idof < rdof; idof++)
1007 : 0 : gradp[icomp][idir] += P(e, mark+idof) * dBdxi[idir][idof];
1008 : : }
1009 : : }
1010 : :
1011 : : // Store the extrema for the gradients
1012 [ + + ]: 10741830 : for (std::size_t c=0; c<ncomp; ++c)
1013 : : {
1014 [ + + ]: 35521500 : for (std::size_t idof = 0; idof < ndof_NodalExtrm; idof++)
1015 : : {
1016 : 26641125 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1017 : 26641125 : auto min_mark = max_mark + 1;
1018 [ + + ]: 26641125 : auto& ex = uNodalExtrm[i->second];
1019 [ + + ][ + + ]: 33458434 : ex[max_mark] = std::max(ex[max_mark], gradc[c][idof]);
1020 : 26641125 : ex[min_mark] = std::min(ex[min_mark], gradc[c][idof]);
1021 : : }
1022 : : }
1023 [ - + ]: 1861455 : for (std::size_t c=0; c<nprim; ++c)
1024 : : {
1025 [ - - ]: 0 : for (std::size_t idof = 0; idof < ndof_NodalExtrm; idof++)
1026 : : {
1027 : 0 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1028 : 0 : auto min_mark = max_mark + 1;
1029 [ - - ]: 0 : auto& ex = pNodalExtrm[i->second];
1030 [ - - ][ - - ]: 0 : ex[max_mark] = std::max(ex[max_mark], gradp[c][idof]);
1031 : 0 : ex[min_mark] = std::min(ex[min_mark], gradp[c][idof]);
1032 : : }
1033 : : }
1034 : : }
1035 : : }
1036 : : }
1037 : 7485 : }
1038 : :
1039 : : void
1040 : 72480 : DG::lim()
1041 : : // *****************************************************************************
1042 : : // Compute limiter function
1043 : : // *****************************************************************************
1044 : : {
1045 : 72480 : auto d = Disc();
1046 : 72480 : const auto rdof = g_inputdeck.get< tag::rdof >();
1047 : 72480 : const auto ncomp = m_u.nprop() / rdof;
1048 : 72480 : const auto nprim = m_p.nprop() / rdof;
1049 : :
1050 : : // Combine own and communicated contributions to nodal extrema
1051 [ + + ]: 2360640 : for (const auto& [gid,g] : m_uNodalExtrmc) {
1052 : 2288160 : auto bid = tk::cref_find( d->Bid(), gid );
1053 [ + + ]: 8291580 : for (ncomp_t c=0; c<ncomp; ++c)
1054 : : {
1055 [ + + ]: 24013680 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
1056 : : {
1057 : 18010260 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1058 : 18010260 : auto min_mark = max_mark + 1;
1059 : 18010260 : m_uNodalExtrm[bid][max_mark] =
1060 [ + + ][ + + ]: 18773425 : std::max(g[max_mark], m_uNodalExtrm[bid][max_mark]);
1061 : 18010260 : m_uNodalExtrm[bid][min_mark] =
1062 : 18010260 : std::min(g[min_mark], m_uNodalExtrm[bid][min_mark]);
1063 : : }
1064 : : }
1065 : : }
1066 [ + + ]: 2360640 : for (const auto& [gid,g] : m_pNodalExtrmc) {
1067 : 2288160 : auto bid = tk::cref_find( d->Bid(), gid );
1068 [ + + ]: 3494370 : for (ncomp_t c=0; c<nprim; ++c)
1069 : : {
1070 [ + + ]: 4824840 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
1071 : : {
1072 : 3618630 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1073 : 3618630 : auto min_mark = max_mark + 1;
1074 : 3618630 : m_pNodalExtrm[bid][max_mark] =
1075 [ - + ][ - + ]: 3618630 : std::max(g[max_mark], m_pNodalExtrm[bid][max_mark]);
1076 : 3618630 : m_pNodalExtrm[bid][min_mark] =
1077 : 3618630 : std::min(g[min_mark], m_pNodalExtrm[bid][min_mark]);
1078 : : }
1079 : : }
1080 : : }
1081 : :
1082 : : // clear gradients receive buffer
1083 : 72480 : tk::destroy(m_uNodalExtrmc);
1084 : 72480 : tk::destroy(m_pNodalExtrmc);
1085 : :
1086 [ + + ]: 72480 : if (rdof > 1) {
1087 : 89280 : g_dgpde[d->MeshId()].limit( d->T(), myGhosts()->m_geoFace, myGhosts()->m_geoElem,
1088 : 44640 : myGhosts()->m_fd, myGhosts()->m_esup, myGhosts()->m_inpoel,
1089 : 44640 : myGhosts()->m_coord, m_ndof, d->Gid(), d->Bid(), m_uNodalExtrm,
1090 : 44640 : m_pNodalExtrm, m_mtInv, m_u, m_p, m_shockmarker );
1091 : :
1092 [ + + ]: 44640 : if (g_inputdeck.get< tag::limsol_projection >())
1093 : 81030 : g_dgpde[d->MeshId()].CPL(m_p, myGhosts()->m_geoElem,
1094 : 40515 : myGhosts()->m_inpoel, myGhosts()->m_coord, m_u,
1095 : 40515 : myGhosts()->m_fd.Esuel().size()/4);
1096 : : }
1097 : :
1098 : : // Send limited solution to neighboring chares
1099 [ + + ]: 72480 : if (myGhosts()->m_sendGhost.empty())
1100 : 4335 : comlim_complete();
1101 : : else
1102 [ + + ]: 636060 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
1103 [ + - ]: 499770 : std::vector< std::size_t > tetid( ghostdata.size() );
1104 [ + - ][ + - ]: 999540 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
[ - - ]
1105 [ + - ]: 499770 : prim( ghostdata.size() );
1106 : : std::vector< std::size_t > ndof;
1107 : : std::size_t j = 0;
1108 [ + + ]: 10020420 : for(const auto& i : ghostdata) {
1109 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4,
1110 : : "Sending limiter ghost data" );
1111 [ + - ]: 9520650 : tetid[j] = i;
1112 [ + - ][ - + ]: 9520650 : u[j] = m_u[i];
1113 [ + - ][ - + ]: 9520650 : prim[j] = m_p[i];
1114 : 9520650 : ++j;
1115 : : }
1116 [ + - ][ + - ]: 999540 : thisProxy[ cid ].comlim( thisIndex, tetid, u, prim );
1117 : : }
1118 : :
1119 : 72480 : ownlim_complete();
1120 : 72480 : }
1121 : :
1122 : : void
1123 : 1770 : DG::refine_ndof()
1124 : : // *****************************************************************************
1125 : : // p-refine all elements that are adjacent to p-refined elements
1126 : : //! \details This function p-refines all the neighbors of an element that has
1127 : : //! been p-refined as a result of an error indicator.
1128 : : // *****************************************************************************
1129 : : {
1130 : 1770 : auto d = Disc();
1131 : 1770 : const auto& inpoel = d->Inpoel();
1132 : 1770 : const auto npoin = d->Coord()[0].size();
1133 : 1770 : const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
1134 : 1770 : std::vector<std::size_t> node_ndof(npoin, 1);
1135 : :
1136 : : // Mark the max ndof for each node and store in node_ndof
1137 [ + + ]: 188540 : for(std::size_t ip=0; ip<npoin; ip++)
1138 : : {
1139 [ + - ]: 186770 : const auto& pesup = tk::cref_find(myGhosts()->m_esup, ip);
1140 [ + + ]: 2641810 : for(auto er : pesup)
1141 [ + + ]: 2537663 : node_ndof[ip] = std::max(m_ndof[er], node_ndof[ip]);
1142 : : }
1143 : :
1144 [ + + ]: 414490 : for(std::size_t e = 0; e < nelem; e++)
1145 : : {
1146 : : // Find if any node of this element has p1/p2 ndofs
1147 : : std::size_t counter_p2(0);
1148 : : std::size_t counter_p1(0);
1149 [ + + ]: 2063600 : for(std::size_t inode = 0; inode < 4; inode++)
1150 : : {
1151 [ + + ]: 1650880 : auto node = inpoel[4*e+inode];
1152 [ + + ]: 1650880 : if(node_ndof[node] == 10)
1153 : 320808 : counter_p2++;
1154 [ + + ]: 1330072 : else if (node_ndof[node] == 4)
1155 : 180140 : counter_p1++;
1156 : : }
1157 : :
1158 : : // If there is at least one node with p1/p2 ndofs, all of the elements
1159 : : // around this node are refined to p1/p2.
1160 [ + + ][ + + ]: 412720 : if(counter_p2 > 0 && m_ndof[e] < 10)
1161 : : {
1162 [ + + ]: 16717 : if(m_ndof[e] == 4)
1163 : 15693 : m_ndof[e] = 10;
1164 [ + + ]: 16717 : if(m_ndof[e] == 1)
1165 : 1024 : m_ndof[e] = 4;
1166 : : }
1167 [ + + ][ + + ]: 396003 : else if(counter_p1 > 0 && m_ndof[e] < 4)
1168 : 13452 : m_ndof[e] = 4;
1169 : : }
1170 : 1770 : }
1171 : :
1172 : 1770 : void DG::smooth_ndof()
1173 : : // *****************************************************************************
1174 : : // Smooth the refined ndof distribution to avoid zigzag refinement
1175 : : // *****************************************************************************
1176 : : {
1177 : 1770 : auto d = Disc();
1178 : 1770 : const auto& inpoel = d->Inpoel();
1179 : 1770 : const auto npoin = d->Coord()[0].size();
1180 : 1770 : const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
1181 : 1770 : std::vector<std::size_t> node_ndof(npoin, 1);
1182 : :
1183 : : // Mark the max ndof for each node and store in node_ndof
1184 [ + + ]: 188540 : for(std::size_t ip=0; ip<npoin; ip++)
1185 : : {
1186 [ + - ]: 186770 : const auto& pesup = tk::cref_find(myGhosts()->m_esup, ip);
1187 [ + + ]: 2641810 : for(auto er : pesup)
1188 [ + + ]: 2543040 : node_ndof[ip] = std::max(m_ndof[er], node_ndof[ip]);
1189 : : }
1190 : :
1191 [ + + ]: 414490 : for(std::size_t e = 0; e < nelem; e++)
1192 : : {
1193 : : // Find if any node of this element has p1/p2 ndofs
1194 : : std::size_t counter_p2(0);
1195 : : std::size_t counter_p1(0);
1196 [ + + ]: 2063600 : for(std::size_t inode = 0; inode < 4; inode++)
1197 : : {
1198 [ + + ]: 1650880 : auto node = inpoel[4*e+inode];
1199 [ + + ]: 1650880 : if(node_ndof[node] == 10)
1200 : 382503 : counter_p2++;
1201 [ + + ]: 1268377 : else if (node_ndof[node] == 4)
1202 : 182564 : counter_p1++;
1203 : : }
1204 : :
1205 : : // If all the nodes in the element are p1/p2, this element is refined to
1206 : : // p1/p2.
1207 [ + + ][ + + ]: 412720 : if(counter_p2 == 4 && m_ndof[e] == 4)
1208 : 1469 : m_ndof[e] = 10;
1209 [ + + ][ + + ]: 411251 : else if(counter_p1 == 4 && m_ndof[e] == 1)
1210 : 1553 : m_ndof[e] = 4;
1211 : : }
1212 : 1770 : }
1213 : :
1214 : : void
1215 : 499770 : DG::comlim( int fromch,
1216 : : const std::vector< std::size_t >& tetid,
1217 : : const std::vector< std::vector< tk::real > >& u,
1218 : : const std::vector< std::vector< tk::real > >& prim )
1219 : : // *****************************************************************************
1220 : : // Receive chare-boundary limiter ghost data from neighboring chares
1221 : : //! \param[in] fromch Sender chare id
1222 : : //! \param[in] tetid Ghost tet ids we receive solution data for
1223 : : //! \param[in] u Limited high-order solution
1224 : : //! \param[in] prim Limited high-order primitive quantities
1225 : : //! \details This function receives contributions to the limited solution from
1226 : : //! fellow chares.
1227 : : // *****************************************************************************
1228 : : {
1229 : : Assert( u.size() == tetid.size(), "Size mismatch in DG::comlim()" );
1230 : : Assert( prim.size() == tetid.size(), "Size mismatch in DG::comlim()" );
1231 : :
1232 : : // Find local-to-ghost tet id map for sender chare
1233 : 499770 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
1234 : :
1235 [ + + ]: 10020420 : for (std::size_t i=0; i<tetid.size(); ++i) {
1236 : 9520650 : auto j = tk::cref_find( n, tetid[i] );
1237 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
1238 : : "Receiving solution non-ghost data" );
1239 : 9520650 : auto b = tk::cref_find( myGhosts()->m_bid, j );
1240 : : Assert( b < m_uc[2].size(), "Indexing out of bounds" );
1241 : : Assert( b < m_pc[2].size(), "Indexing out of bounds" );
1242 : 9520650 : m_uc[2][b] = u[i];
1243 : 9520650 : m_pc[2][b] = prim[i];
1244 : : }
1245 : :
1246 : : // if we have received all solution ghost contributions from neighboring
1247 : : // chares (chares we communicate along chare-boundary faces with), and
1248 : : // contributed our solution to these neighbors, proceed to limiting
1249 [ + + ]: 499770 : if (++m_nlim == myGhosts()->m_sendGhost.size()) {
1250 : 68145 : m_nlim = 0;
1251 : 68145 : comlim_complete();
1252 : : }
1253 : 499770 : }
1254 : :
1255 : : void
1256 : 72480 : DG::dt()
1257 : : // *****************************************************************************
1258 : : // Compute time step size
1259 : : // *****************************************************************************
1260 : : {
1261 : 72480 : auto d = Disc();
1262 : :
1263 : : // Combine own and communicated contributions of limited solution and degrees
1264 : : // of freedom in cells (if p-adaptive)
1265 [ + + ]: 9665610 : for (const auto& b : myGhosts()->m_bid) {
1266 : : Assert( m_uc[2][b.second].size() == m_u.nprop(), "ncomp size mismatch" );
1267 : : Assert( m_pc[2][b.second].size() == m_p.nprop(), "ncomp size mismatch" );
1268 [ + + ]: 149542125 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
1269 : 140021475 : m_u(b.first,c) = m_uc[2][b.second][c];
1270 : : }
1271 [ + + ]: 30332715 : for (std::size_t c=0; c<m_p.nprop(); ++c) {
1272 : 20812065 : m_p(b.first,c) = m_pc[2][b.second][c];
1273 : : }
1274 : : }
1275 : :
1276 : 72480 : auto mindt = std::numeric_limits< tk::real >::max();
1277 : :
1278 [ + + ]: 72480 : if (m_stage == 0)
1279 : : {
1280 [ + + ]: 24160 : auto const_dt = g_inputdeck.get< tag::dt >();
1281 : : auto eps = std::numeric_limits< tk::real >::epsilon();
1282 : :
1283 : : // use constant dt if configured
1284 [ + + ]: 24160 : if (std::abs(const_dt) > eps) {
1285 : :
1286 : 21755 : mindt = const_dt;
1287 : :
1288 : : } else { // compute dt based on CFL
1289 : :
1290 : : // find the minimum dt across all PDEs integrated
1291 : : auto eqdt =
1292 : 4810 : g_dgpde[d->MeshId()].dt( myGhosts()->m_coord, myGhosts()->m_inpoel,
1293 : 2405 : myGhosts()->m_fd,
1294 : 2405 : myGhosts()->m_geoFace, myGhosts()->m_geoElem, m_ndof, m_u, m_p,
1295 : 2405 : myGhosts()->m_fd.Esuel().size()/4 );
1296 [ + - ]: 2405 : if (eqdt < mindt) mindt = eqdt;
1297 : :
1298 : 2405 : mindt *= g_inputdeck.get< tag::cfl >();
1299 : : }
1300 : : }
1301 : : else
1302 : : {
1303 : 48320 : mindt = d->Dt();
1304 : : }
1305 : :
1306 : : // Resize the buffer vector of nodal extrema
1307 : 72480 : resizeNodalExtremac();
1308 : :
1309 : : // Contribute to minimum dt across all chares then advance to next step
1310 [ + - ]: 72480 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
1311 : 72480 : CkCallback(CkReductionTarget(DG,solve), thisProxy) );
1312 : 72480 : }
1313 : :
1314 : : void
1315 : 72480 : DG::solve( tk::real newdt )
1316 : : // *****************************************************************************
1317 : : // Compute right-hand side of discrete transport equations
1318 : : //! \param[in] newdt Size of this new time step
1319 : : // *****************************************************************************
1320 : : {
1321 : : // Enable SDAG wait for building the solution vector during the next stage
1322 [ + - ]: 72480 : thisProxy[ thisIndex ].wait4sol();
1323 [ + - ]: 72480 : thisProxy[ thisIndex ].wait4refine();
1324 [ + - ]: 72480 : thisProxy[ thisIndex ].wait4smooth();
1325 [ + - ]: 72480 : thisProxy[ thisIndex ].wait4reco();
1326 [ + - ]: 72480 : thisProxy[ thisIndex ].wait4nodalExtrema();
1327 [ + - ]: 72480 : thisProxy[ thisIndex ].wait4lim();
1328 [ + - ]: 72480 : thisProxy[ thisIndex ].wait4nod();
1329 : :
1330 : 72480 : auto d = Disc();
1331 : 72480 : const auto rdof = g_inputdeck.get< tag::rdof >();
1332 : 72480 : const auto ndof = g_inputdeck.get< tag::ndof >();
1333 : 72480 : const auto neq = m_u.nprop()/rdof;
1334 : :
1335 : : // Set new time step size
1336 [ + + ]: 72480 : if (m_stage == 0) d->setdt( newdt );
1337 : :
1338 : 72480 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
1339 [ + + ][ + + ]: 72480 : if (pref && m_stage == 0)
1340 : : {
1341 : : // When the element are coarsened, high order terms should be zero
1342 [ + + ]: 809600 : for(std::size_t e = 0; e < myGhosts()->m_nunk; e++)
1343 : : {
1344 : 807830 : const auto ncomp= m_u.nprop()/rdof;
1345 [ + + ]: 807830 : if(m_ndof[e] == 1)
1346 : : {
1347 [ + + ]: 3427608 : for (std::size_t c=0; c<ncomp; ++c)
1348 : : {
1349 : 2856340 : auto mark = c*rdof;
1350 : 2856340 : m_u(e, mark+1) = 0.0;
1351 : 2856340 : m_u(e, mark+2) = 0.0;
1352 : 2856340 : m_u(e, mark+3) = 0.0;
1353 : : }
1354 [ + + ]: 236562 : } else if(m_ndof[e] == 4)
1355 : : {
1356 [ + + ]: 508752 : for (std::size_t c=0; c<ncomp; ++c)
1357 : : {
1358 : 423960 : auto mark = c*ndof;
1359 : 423960 : m_u(e, mark+4) = 0.0;
1360 : 423960 : m_u(e, mark+5) = 0.0;
1361 : 423960 : m_u(e, mark+6) = 0.0;
1362 : 423960 : m_u(e, mark+7) = 0.0;
1363 : 423960 : m_u(e, mark+8) = 0.0;
1364 : 423960 : m_u(e, mark+9) = 0.0;
1365 : : }
1366 : : }
1367 : : }
1368 : : }
1369 : :
1370 : : // Update Un
1371 [ + + ]: 72480 : if (m_stage == 0) m_un = m_u;
1372 : :
1373 : : // physical time at time-stage for computing exact source terms
1374 : 72480 : tk::real physT(d->T());
1375 [ + + ]: 72480 : if (m_stage == 1) {
1376 : 24160 : physT += d->Dt();
1377 : : }
1378 [ + + ]: 48320 : else if (m_stage == 2) {
1379 : 24160 : physT += 0.5*d->Dt();
1380 : : }
1381 : :
1382 : 144960 : g_dgpde[d->MeshId()].rhs( physT, myGhosts()->m_geoFace, myGhosts()->m_geoElem,
1383 : 72480 : myGhosts()->m_fd, myGhosts()->m_inpoel, m_boxelems, myGhosts()->m_coord,
1384 : 72480 : m_u, m_p, m_ndof, d->Dt(), m_rhs );
1385 : :
1386 : : // Explicit time-stepping using RK3 to discretize time-derivative
1387 [ + + ]: 28217460 : for(std::size_t e=0; e<myGhosts()->m_nunk; ++e)
1388 [ + + ]: 118000080 : for(std::size_t c=0; c<neq; ++c)
1389 : : {
1390 [ + + ]: 383278230 : for (std::size_t k=0; k<m_numEqDof[c]; ++k)
1391 : : {
1392 : 293423130 : auto rmark = c*rdof+k;
1393 : 293423130 : auto mark = c*ndof+k;
1394 [ + + ]: 293423130 : m_u(e, rmark) = rkcoef[0][m_stage] * m_un(e, rmark)
1395 : 293423130 : + rkcoef[1][m_stage] * ( m_u(e, rmark)
1396 : 293423130 : + d->Dt() * m_rhs(e, mark)/m_lhs(e, mark) );
1397 [ + + ]: 293423130 : if(fabs(m_u(e, rmark)) < 1e-16)
1398 : 135731391 : m_u(e, rmark) = 0;
1399 : : }
1400 : : // zero out unused/reconstructed dofs of equations using reduced dofs
1401 : : // (see DGMultiMat::numEquationDofs())
1402 [ + + ]: 89855100 : if (m_numEqDof[c] < rdof) {
1403 [ + + ]: 77428860 : for (std::size_t k=m_numEqDof[c]; k<rdof; ++k)
1404 : : {
1405 : 58071645 : auto rmark = c*rdof+k;
1406 : 58071645 : m_u(e, rmark) = 0.0;
1407 : : }
1408 : : }
1409 : : }
1410 : :
1411 : : // Update primitives based on the evolved solution
1412 : 72480 : g_dgpde[d->MeshId()].updateInterfaceCells( m_u,
1413 : 72480 : myGhosts()->m_fd.Esuel().size()/4, m_ndof );
1414 : 72480 : g_dgpde[d->MeshId()].updatePrimitives( m_u, m_lhs, myGhosts()->m_geoElem, m_p,
1415 : 72480 : myGhosts()->m_fd.Esuel().size()/4 );
1416 [ + - ]: 72480 : if (!g_inputdeck.get< tag::accuracy_test >()) {
1417 : 72480 : g_dgpde[d->MeshId()].cleanTraceMaterial( physT, myGhosts()->m_geoElem, m_u,
1418 : 72480 : m_p, myGhosts()->m_fd.Esuel().size()/4 );
1419 : : }
1420 : :
1421 [ + + ]: 72480 : if (m_stage < 2) {
1422 : :
1423 : : // continue with next time step stage
1424 : 48320 : stage();
1425 : :
1426 : : } else {
1427 : :
1428 : : // Increase number of iterations and physical time
1429 : 24160 : d->next();
1430 : :
1431 : : // Compute diagnostics, e.g., residuals
1432 : 24160 : auto diag_computed = m_diag.compute( *d,
1433 : 24160 : m_u.nunk()-myGhosts()->m_fd.Esuel().size()/4, myGhosts()->m_geoElem,
1434 : 24160 : m_ndof, m_u, m_un );
1435 : :
1436 : : // Continue to mesh refinement (if configured)
1437 [ + + ][ + - ]: 35200 : if (!diag_computed) refine( std::vector< tk::real >( m_u.nprop(), 0.0 ) );
[ + - ]
1438 : :
1439 : : }
1440 : 72480 : }
1441 : :
1442 : : void
1443 : 24160 : DG::refine( [[maybe_unused]] const std::vector< tk::real >& l2res )
1444 : : // *****************************************************************************
1445 : : // Optionally refine/derefine mesh
1446 : : //! \param[in] l2res L2-norms of the residual for each scalar component
1447 : : //! computed across the whole problem
1448 : : // *****************************************************************************
1449 : : {
1450 : 24160 : auto d = Disc();
1451 : :
1452 : 24160 : auto dtref = g_inputdeck.get< tag::amr, tag::dtref >();
1453 : 24160 : auto dtfreq = g_inputdeck.get< tag::amr, tag::dtfreq >();
1454 : :
1455 : : // if t>0 refinement enabled and we hit the dtref frequency
1456 [ - + ][ - - ]: 24160 : if (dtref && !(d->It() % dtfreq)) { // refine
1457 : :
1458 : 0 : d->startvol();
1459 [ - - ][ - - ]: 0 : d->Ref()->dtref( myGhosts()->m_fd.Bface(), {},
[ - - ][ - - ]
1460 : 0 : tk::remap(myGhosts()->m_fd.Triinpoel(),d->Gid()) );
1461 : 0 : d->refined() = 1;
1462 : :
1463 : : } else { // do not refine
1464 : :
1465 : 24160 : d->refined() = 0;
1466 : 24160 : stage();
1467 : :
1468 : : }
1469 : 24160 : }
1470 : :
1471 : : void
1472 : 0 : DG::resizePostAMR(
1473 : : const std::vector< std::size_t >& /*ginpoel*/,
1474 : : const tk::UnsMesh::Chunk& chunk,
1475 : : const tk::UnsMesh::Coords& coord,
1476 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /*addedNodes*/,
1477 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
1478 : : const std::set< std::size_t >& removedNodes,
1479 : : const std::unordered_map< std::size_t, std::size_t >& amrNodeMap,
1480 : : const tk::NodeCommMap& nodeCommMap,
1481 : : const std::map< int, std::vector< std::size_t > >& bface,
1482 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
1483 : : const std::vector< std::size_t >& triinpoel,
1484 : : const std::unordered_map< std::size_t, std::set< std::size_t > >& elemblockid )
1485 : : // *****************************************************************************
1486 : : // Receive new mesh from Refiner
1487 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
1488 : : //! \param[in] coord New mesh node coordinates
1489 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1490 : : //! \param[in] removedNodes Newly removed mesh node local ids
1491 : : //! \param[in] amrNodeMap Node id map after amr (local ids)
1492 : : //! \param[in] nodeCommMap New node communication map
1493 : : //! \param[in] bface Boundary-faces mapped to side set ids
1494 : : //! \param[in] triinpoel Boundary-face connectivity
1495 : : //! \param[in] elemblockid Local tet ids associated with mesh block ids
1496 : : // *****************************************************************************
1497 : : {
1498 : 0 : auto d = Disc();
1499 : :
1500 : : // Set flag that indicates that we are during time stepping
1501 : 0 : m_initial = 0;
1502 : 0 : myGhosts()->m_initial = 0;
1503 : :
1504 : : // Zero field output iteration count between two mesh refinement steps
1505 : 0 : d->Itf() = 0;
1506 : :
1507 : : // Increase number of iterations with mesh refinement
1508 : 0 : ++d->Itr();
1509 : :
1510 : : // Save old number of elements
1511 : 0 : [[maybe_unused]] auto old_nelem = myGhosts()->m_inpoel.size()/4;
1512 : :
1513 : : // Resize mesh data structures
1514 : 0 : d->resizePostAMR( chunk, coord, amrNodeMap, nodeCommMap, removedNodes,
1515 : : elemblockid );
1516 : :
1517 : : // Update state
1518 : 0 : myGhosts()->m_inpoel = d->Inpoel();
1519 : 0 : myGhosts()->m_coord = d->Coord();
1520 : 0 : auto nelem = myGhosts()->m_inpoel.size()/4;
1521 : : m_p.resize( nelem );
1522 : : m_u.resize( nelem );
1523 : : m_un.resize( nelem );
1524 : : m_lhs.resize( nelem );
1525 : : m_rhs.resize( nelem );
1526 [ - - ][ - - ]: 0 : m_uNodalExtrm.resize( Disc()->Bid().size(), std::vector<tk::real>( 2*
1527 [ - - ]: 0 : m_ndof_NodalExtrm*g_inputdeck.get< tag::ncomp >() ) );
1528 [ - - ][ - - ]: 0 : m_pNodalExtrm.resize( Disc()->Bid().size(), std::vector<tk::real>( 2*
1529 [ - - ]: 0 : m_ndof_NodalExtrm*m_p.nprop()/g_inputdeck.get< tag::rdof >()));
1530 : :
1531 : : // Resize the buffer vector of nodal extrema
1532 : 0 : resizeNodalExtremac();
1533 : :
1534 [ - - ][ - - ]: 0 : myGhosts()->m_fd = FaceData( myGhosts()->m_inpoel, bface,
[ - - ][ - - ]
[ - - ]
1535 : 0 : tk::remap(triinpoel,d->Lid()) );
1536 : :
1537 [ - - ]: 0 : myGhosts()->m_geoFace =
1538 : 0 : tk::Fields( tk::genGeoFaceTri( myGhosts()->m_fd.Nipfac(),
1539 : 0 : myGhosts()->m_fd.Inpofa(), coord ) );
1540 [ - - ]: 0 : myGhosts()->m_geoElem = tk::Fields( tk::genGeoElemTet( myGhosts()->m_inpoel,
1541 : 0 : coord ) );
1542 : :
1543 : 0 : myGhosts()->m_nfac = myGhosts()->m_fd.Inpofa().size()/3;
1544 : 0 : myGhosts()->m_nunk = nelem;
1545 : 0 : m_npoin = coord[0].size();
1546 : 0 : myGhosts()->m_bndFace.clear();
1547 : 0 : myGhosts()->m_exptGhost.clear();
1548 : 0 : myGhosts()->m_sendGhost.clear();
1549 : 0 : myGhosts()->m_ghost.clear();
1550 : 0 : myGhosts()->m_esup.clear();
1551 : :
1552 : : // Update solution on new mesh, P0 (cell center value) only for now
1553 : : m_un = m_u;
1554 : : auto pn = m_p;
1555 : 0 : auto unprop = m_u.nprop();
1556 : : auto pnprop = m_p.nprop();
1557 [ - - ]: 0 : for (const auto& [child,parent] : addedTets) {
1558 : : Assert( child < nelem, "Indexing out of new solution vector" );
1559 : : Assert( parent < old_nelem, "Indexing out of old solution vector" );
1560 [ - - ]: 0 : for (std::size_t i=0; i<unprop; ++i) m_u(child,i) = m_un(parent,i);
1561 [ - - ]: 0 : for (std::size_t i=0; i<pnprop; ++i) m_p(child,i) = pn(parent,i);
1562 : : }
1563 : : m_un = m_u;
1564 : :
1565 : : // Resize communication buffers
1566 [ - - ][ - - ]: 0 : m_ghosts[thisIndex].resizeComm();
[ - - ][ - - ]
1567 : 0 : }
1568 : :
1569 : : bool
1570 : 21801 : DG::fieldOutput() const
1571 : : // *****************************************************************************
1572 : : // Decide wether to output field data
1573 : : //! \return True if field data is output in this step
1574 : : // *****************************************************************************
1575 : : {
1576 : 21801 : auto d = Disc();
1577 : :
1578 : : // Output field data
1579 [ + + ][ + - ]: 21801 : return d->fielditer() or d->fieldtime() or d->fieldrange() or d->finished();
[ + - ][ + + ]
1580 : : }
1581 : :
1582 : : bool
1583 : 2176 : DG::refinedOutput() const
1584 : : // *****************************************************************************
1585 : : // Decide if we write field output using a refined mesh
1586 : : //! \return True if field output will use a refined mesh
1587 : : // *****************************************************************************
1588 : : {
1589 [ + + ]: 2176 : return g_inputdeck.get< tag::field_output, tag::refined >() &&
1590 [ - + ]: 33 : g_inputdeck.get< tag::scheme >() != ctr::SchemeType::DG;
1591 : : }
1592 : :
1593 : : void
1594 : 2176 : DG::writeFields(
1595 : : CkCallback c,
1596 : : const std::unordered_map< std::size_t, std::size_t >& addedTets )
1597 : : // *****************************************************************************
1598 : : // Output mesh field data
1599 : : //! \param[in] c Function to continue with after the write
1600 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1601 : : // *****************************************************************************
1602 : : {
1603 : 2176 : auto d = Disc();
1604 : :
1605 : : const auto& inpoel = std::get< 0 >( m_outmesh.chunk );
1606 : 4352 : auto esup = tk::genEsup( inpoel, 4 );
1607 : 2176 : auto nelem = inpoel.size() / 4;
1608 : :
1609 : : // Combine own and communicated contributions and finish averaging of node
1610 : : // field output in chare boundary nodes
1611 : : const auto& lid = std::get< 2 >( m_outmesh.chunk );
1612 [ + + ]: 87374 : for (const auto& [g,f] : m_uNodefieldsc) {
1613 : : Assert( m_uNodefields.nprop() == f.first.size(), "Size mismatch" );
1614 : 85198 : auto p = tk::cref_find( lid, g );
1615 [ + + ]: 347454 : for (std::size_t i=0; i<f.first.size(); ++i) {
1616 : 262256 : m_uNodefields(p,i) += f.first[i];
1617 : 262256 : m_uNodefields(p,i) /= static_cast< tk::real >(
1618 : 262256 : esup.second[p+1] - esup.second[p] + f.second );
1619 : : }
1620 : : }
1621 : 2176 : tk::destroy( m_uNodefieldsc );
1622 [ + + ]: 87374 : for (const auto& [g,f] : m_pNodefieldsc) {
1623 : : Assert( m_pNodefields.nprop() == f.first.size(), "Size mismatch" );
1624 : 85198 : auto p = tk::cref_find( lid, g );
1625 [ + + ]: 123582 : for (std::size_t i=0; i<f.first.size(); ++i) {
1626 : 38384 : m_pNodefields(p,i) += f.first[i];
1627 : 38384 : m_pNodefields(p,i) /= static_cast< tk::real >(
1628 : 38384 : esup.second[p+1] - esup.second[p] + f.second );
1629 : : }
1630 : : }
1631 : 2176 : tk::destroy( m_pNodefieldsc );
1632 : :
1633 : : // Lambda to decide if a node (global id) is on a chare boundary of the field
1634 : : // output mesh. p - global node id, return true if node is on the chare
1635 : : // boundary.
1636 : : auto chbnd = [ this ]( std::size_t p ) {
1637 : : return
1638 : : std::any_of( m_outmesh.nodeCommMap.cbegin(), m_outmesh.nodeCommMap.cend(),
1639 : : [&](const auto& s) { return s.second.find(p) != s.second.cend(); } );
1640 : : };
1641 : :
1642 : : // Finish computing node field output averages in internal nodes
1643 : 2176 : auto npoin = m_outmesh.coord[0].size();
1644 : : auto& gid = std::get< 1 >( m_outmesh.chunk );
1645 [ + + ]: 416390 : for (std::size_t p=0; p<npoin; ++p) {
1646 [ + + ]: 414214 : if (!chbnd(gid[p])) {
1647 : 329016 : auto n = static_cast< tk::real >( esup.second[p+1] - esup.second[p] );
1648 [ + + ]: 1320705 : for (std::size_t i=0; i<m_uNodefields.nprop(); ++i)
1649 : 991689 : m_uNodefields(p,i) /= n;
1650 [ + + ]: 458108 : for (std::size_t i=0; i<m_pNodefields.nprop(); ++i)
1651 : 129092 : m_pNodefields(p,i) /= n;
1652 : : }
1653 : : }
1654 : :
1655 : : // Collect field output from numerical solution requested by user
1656 : 2176 : auto elemfields = numericFieldOutput( m_uElemfields, tk::Centering::ELEM,
1657 [ + - ][ + - ]: 4352 : g_dgpde[Disc()->MeshId()].OutVarFn(), m_pElemfields );
[ + - ]
1658 : 2176 : auto nodefields = numericFieldOutput( m_uNodefields, tk::Centering::NODE,
1659 [ + - ][ + - ]: 4352 : g_dgpde[Disc()->MeshId()].OutVarFn(), m_pNodefields );
[ + - ]
1660 : :
1661 : : // Collect field output from analytical solutions (if exist)
1662 : 2176 : const auto& coord = m_outmesh.coord;
1663 [ + - ]: 2176 : auto geoElem = tk::genGeoElemTet( inpoel, coord );
1664 [ + - ]: 2176 : auto t = Disc()->T();
1665 [ + - ]: 2176 : analyticFieldOutput( g_dgpde[d->MeshId()], tk::Centering::ELEM,
1666 [ + - ][ + - ]: 8704 : geoElem.extract_comp(1), geoElem.extract_comp(2), geoElem.extract_comp(3),
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
1667 : : t, elemfields );
1668 [ + - ]: 2176 : analyticFieldOutput( g_dgpde[d->MeshId()], tk::Centering::NODE, coord[0],
1669 : : coord[1], coord[2], t, nodefields );
1670 : :
1671 : : // Add adaptive indicator array to element-centered field output
1672 [ + + ]: 2176 : if (g_inputdeck.get< tag::pref, tag::pref >()) {
1673 [ + - ]: 402 : std::vector< tk::real > ndof( begin(m_ndof), end(m_ndof) );
1674 [ + - ]: 402 : ndof.resize( nelem );
1675 [ + + ]: 87834 : for (const auto& [child,parent] : addedTets)
1676 : 87432 : ndof[child] = static_cast< tk::real >( m_ndof[parent] );
1677 [ + - ]: 402 : elemfields.push_back( ndof );
1678 : : }
1679 : :
1680 : : // Add shock detection marker array to element-centered field output
1681 [ + - ][ - - ]: 2176 : std::vector< tk::real > shockmarker( begin(m_shockmarker), end(m_shockmarker) );
1682 : : // Here m_shockmarker has a size of m_u.nunk() which is the number of the
1683 : : // elements within this partition (nelem) plus the ghost partition cells. In
1684 : : // terms of output purpose, we only need the solution data within this
1685 : : // partition. Therefore, resizing it to nelem removes the extra partition
1686 : : // boundary allocations in the shockmarker vector. Since the code assumes that
1687 : : // the boundary elements are on the top, the resize operation keeps the lower
1688 : : // portion.
1689 [ + - ]: 2176 : shockmarker.resize( nelem );
1690 [ + + ]: 159688 : for (const auto& [child,parent] : addedTets)
1691 : 157512 : shockmarker[child] = static_cast< tk::real >(m_shockmarker[parent]);
1692 [ + - ]: 2176 : elemfields.push_back( shockmarker );
1693 : :
1694 : : // Add inverse deformation gradient tensor to element-centered field output
1695 [ + - ][ + - ]: 2176 : auto defgrad = g_dgpde[d->MeshId()].cellAvgDeformGrad(m_u,
1696 [ + - ][ + - ]: 2176 : myGhosts()->m_fd.Esuel().size()/4);
1697 [ + + ]: 2176 : if (!defgrad[0].empty()) {
1698 [ - + ]: 14 : for (const auto& [child,parent] : addedTets)
1699 [ - - ]: 0 : for (auto& gij : defgrad)
1700 : 0 : gij[child] = static_cast< tk::real >(gij[parent]);
1701 [ + + ]: 140 : for (const auto& gij : defgrad)
1702 [ + - ]: 126 : elemfields.push_back(gij);
1703 : : }
1704 : :
1705 : : // Query fields names requested by user
1706 [ + - ]: 4352 : auto elemfieldnames = numericFieldNames( tk::Centering::ELEM );
1707 [ + - ]: 2176 : auto nodefieldnames = numericFieldNames( tk::Centering::NODE );
1708 : :
1709 : : // Collect field output names for analytical solutions
1710 [ + - ]: 2176 : analyticFieldNames( g_dgpde[d->MeshId()], tk::Centering::ELEM, elemfieldnames );
1711 [ + - ]: 2176 : analyticFieldNames( g_dgpde[d->MeshId()], tk::Centering::NODE, nodefieldnames );
1712 : :
1713 [ + + ]: 2176 : if (g_inputdeck.get< tag::pref, tag::pref >())
1714 [ + - ]: 804 : elemfieldnames.push_back( "NDOF" );
1715 : :
1716 [ + - ][ + + ]: 4352 : elemfieldnames.push_back( "shock_marker" );
1717 : :
1718 [ + + ]: 2176 : if (!defgrad[0].empty()) {
1719 [ + + ]: 56 : for (std::size_t i=1; i<=3; ++i)
1720 [ + + ]: 168 : for (std::size_t j=1; j<=3; ++j)
1721 [ + - ][ + - ]: 252 : elemfieldnames.push_back("g"+std::to_string(i)+std::to_string(j));
[ - + ][ - + ]
[ - + ][ - - ]
[ - - ][ - - ]
1722 : : }
1723 : :
1724 : : Assert( elemfieldnames.size() == elemfields.size(), "Size mismatch" );
1725 : : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
1726 : :
1727 : : // Output chare mesh and fields metadata to file
1728 : 2176 : const auto& triinpoel = m_outmesh.triinpoel;
1729 [ + - ][ + - ]: 6528 : d->write( inpoel, m_outmesh.coord, m_outmesh.bface, {},
[ + + ][ - - ]
1730 [ + - ]: 4352 : tk::remap( triinpoel, lid ), elemfieldnames, nodefieldnames,
1731 : : {}, {}, elemfields, nodefields, {}, {}, c );
1732 : 2176 : }
1733 : :
1734 : : void
1735 : 13330 : DG::comnodeout( const std::vector< std::size_t >& gid,
1736 : : const std::vector< std::size_t >& nesup,
1737 : : const std::vector< std::vector< tk::real > >& Lu,
1738 : : const std::vector< std::vector< tk::real > >& Lp )
1739 : : // *****************************************************************************
1740 : : // Receive chare-boundary nodal solution (for field output) contributions from
1741 : : // neighboring chares
1742 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
1743 : : //! \param[in] nesup Number of elements surrounding points
1744 : : //! \param[in] Lu Partial contributions of solution nodal fields to
1745 : : //! chare-boundary nodes
1746 : : //! \param[in] Lp Partial contributions of primitive quantity nodal fields to
1747 : : //! chare-boundary nodes
1748 : : // *****************************************************************************
1749 : : {
1750 : : Assert( gid.size() == nesup.size(), "Size mismatch" );
1751 : : Assert(Lu.size() == m_uNodefields.nprop(), "Fields size mismatch");
1752 : : Assert(Lp.size() == m_pNodefields.nprop(), "Fields size mismatch");
1753 [ + + ]: 62724 : for (std::size_t f=0; f<Lu.size(); ++f)
1754 : : Assert( gid.size() == Lu[f].size(), "Size mismatch" );
1755 [ + + ]: 21242 : for (std::size_t f=0; f<Lp.size(); ++f)
1756 : : Assert( gid.size() == Lp[f].size(), "Size mismatch" );
1757 : :
1758 [ + + ]: 141372 : for (std::size_t i=0; i<gid.size(); ++i) {
1759 : : auto& nfu = m_uNodefieldsc[ gid[i] ];
1760 : 128042 : nfu.first.resize( Lu.size() );
1761 [ + + ]: 587774 : for (std::size_t f=0; f<Lu.size(); ++f) nfu.first[f] += Lu[f][i];
1762 : 128042 : nfu.second += nesup[i];
1763 : 128042 : auto& nfp = m_pNodefieldsc[ gid[i] ];
1764 : 128042 : nfp.first.resize( Lp.size() );
1765 [ + + ]: 209066 : for (std::size_t f=0; f<Lp.size(); ++f) nfp.first[f] += Lp[f][i];
1766 : 128042 : nfp.second += nesup[i];
1767 : : }
1768 : :
1769 : : // When we have heard from all chares we communicate with, this chare is done
1770 [ + + ]: 13330 : if (++m_nnod == Disc()->NodeCommMap().size()) {
1771 : 1996 : m_nnod = 0;
1772 : 1996 : comnodeout_complete();
1773 : : }
1774 : 13330 : }
1775 : :
1776 : : void
1777 : 72480 : DG::stage()
1778 : : // *****************************************************************************
1779 : : // Evaluate whether to continue with next time step stage
1780 : : // *****************************************************************************
1781 : : {
1782 : : // Increment Runge-Kutta stage counter
1783 : 72480 : ++m_stage;
1784 : :
1785 : : // if not all Runge-Kutta stages complete, continue to next time stage,
1786 : : // otherwise prepare for nodal field output
1787 [ + + ]: 72480 : if (m_stage < 3)
1788 : 48320 : next();
1789 : : else
1790 [ + - ][ + - ]: 72480 : startFieldOutput( CkCallback(CkIndex_DG::step(), thisProxy[thisIndex]) );
[ - + ][ - - ]
1791 : 72480 : }
1792 : :
1793 : : void
1794 : 23307 : DG::evalLB( int nrestart )
1795 : : // *****************************************************************************
1796 : : // Evaluate whether to do load balancing
1797 : : //! \param[in] nrestart Number of times restarted
1798 : : // *****************************************************************************
1799 : : {
1800 : 23307 : auto d = Disc();
1801 : :
1802 : : // Detect if just returned from a checkpoint and if so, zero timers
1803 : 23307 : d->restarted( nrestart );
1804 : :
1805 : 23307 : const auto lbfreq = g_inputdeck.get< tag::cmd, tag::lbfreq >();
1806 : 23307 : const auto nonblocking = g_inputdeck.get< tag::cmd, tag::nonblocking >();
1807 : :
1808 : : // Load balancing if user frequency is reached or after the second time-step
1809 [ + + ][ + + ]: 23307 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
1810 : :
1811 : 20187 : AtSync();
1812 [ - + ]: 20187 : if (nonblocking) next();
1813 : :
1814 : : } else {
1815 : :
1816 : 3120 : next();
1817 : :
1818 : : }
1819 : 23307 : }
1820 : :
1821 : : void
1822 : 23307 : DG::evalRestart()
1823 : : // *****************************************************************************
1824 : : // Evaluate whether to save checkpoint/restart
1825 : : // *****************************************************************************
1826 : : {
1827 : 23307 : auto d = Disc();
1828 : :
1829 : 23307 : const auto rsfreq = g_inputdeck.get< tag::cmd, tag::rsfreq >();
1830 : 23307 : const auto benchmark = g_inputdeck.get< tag::cmd, tag::benchmark >();
1831 : :
1832 [ + + ][ - + ]: 23307 : if (not benchmark and not (d->It() % rsfreq)) {
1833 : :
1834 : 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
1835 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
1836 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
[ - - ]
1837 : :
1838 : : } else {
1839 : :
1840 : 23307 : evalLB( /* nrestart = */ -1 );
1841 : :
1842 : : }
1843 : 23307 : }
1844 : :
1845 : : void
1846 : 24160 : DG::step()
1847 : : // *****************************************************************************
1848 : : // Evaluate wether to continue with next time step
1849 : : // *****************************************************************************
1850 : : {
1851 : 24160 : auto d = Disc();
1852 : :
1853 : : // Output time history
1854 [ + - ][ + - ]: 24160 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ]
1855 : 0 : std::vector< std::vector< tk::real > > hist;
1856 [ - - ][ - - ]: 0 : auto h = g_dgpde[d->MeshId()].histOutput( d->Hist(), myGhosts()->m_inpoel,
1857 [ - - ][ - - ]: 0 : myGhosts()->m_coord, m_u, m_p );
1858 [ - - ]: 0 : hist.insert( end(hist), begin(h), end(h) );
1859 [ - - ]: 0 : d->history( std::move(hist) );
1860 : : }
1861 : :
1862 : : // Free memory storing output mesh
1863 : 24160 : m_outmesh.destroy();
1864 : :
1865 : : // Output one-liner status report to screen
1866 : 24160 : d->status();
1867 : : // Reset Runge-Kutta stage counter
1868 : 24160 : m_stage = 0;
1869 : :
1870 : 24160 : const auto term = g_inputdeck.get< tag::term >();
1871 : 24160 : const auto nstep = g_inputdeck.get< tag::nstep >();
1872 : : const auto eps = std::numeric_limits< tk::real >::epsilon();
1873 : :
1874 : : // If neither max iterations nor max time reached, continue, otherwise finish
1875 [ + - ][ + + ]: 24160 : if (std::fabs(d->T()-term) > eps && d->It() < nstep) {
1876 : :
1877 : 23307 : evalRestart();
1878 : :
1879 : : } else {
1880 : :
1881 : 853 : auto meshid = d->MeshId();
1882 [ + - ]: 1706 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1883 : 1706 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
1884 : :
1885 : : }
1886 : 24160 : }
1887 : :
1888 : : #include "NoWarning/dg.def.h"
|
