Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/FV.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 FV advances a system of PDEs with the finite volume scheme
9 : : \details FV advances a system of partial differential equations (PDEs) using
10 : : the finite volume (FV) (on tetrahedron elements).
11 : : \see The documentation in FV.h.
12 : : */
13 : : // *****************************************************************************
14 : :
15 : : #include <algorithm>
16 : : #include <numeric>
17 : : #include <sstream>
18 : :
19 : : #include "FV.hpp"
20 : : #include "Discretization.hpp"
21 : : #include "FVPDE.hpp"
22 : : #include "DiagReducer.hpp"
23 : : #include "DerivedData.hpp"
24 : : #include "ElemDiagnostics.hpp"
25 : : #include "Inciter/InputDeck/InputDeck.hpp"
26 : : #include "Refiner.hpp"
27 : : #include "Limiter.hpp"
28 : : #include "PrefIndicator.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< FVPDE > g_fvpde;
40 : :
41 : : } // inciter::
42 : :
43 : : extern tk::CProxy_ChareStateCollector stateProxy;
44 : :
45 : : using inciter::FV;
46 : :
47 : 299 : FV::FV( const CProxy_Discretization& disc,
48 : : const CProxy_Ghosts& ghostsproxy,
49 : : const std::map< int, std::vector< std::size_t > >& bface,
50 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
51 : 299 : const std::vector< std::size_t >& triinpoel ) :
52 : : m_disc( disc ),
53 : : m_ghosts( ghostsproxy ),
54 : : m_nsol( 0 ),
55 : : m_ninitsol( 0 ),
56 : : m_nlim( 0 ),
57 : : m_nnod( 0 ),
58 [ + - ]: 299 : m_u( Disc()->Inpoel().size()/4,
59 : 299 : g_inputdeck.get< tag::rdof >()*
60 : 299 : g_inputdeck.get< tag::ncomp >() ),
61 : : m_un( m_u.nunk(), m_u.nprop() ),
62 : 299 : m_p( m_u.nunk(), g_inputdeck.get< tag::rdof >()*
63 [ + - ][ + - ]: 299 : g_fvpde[Disc()->MeshId()].nprim() ),
64 : : m_rhs( m_u.nunk(),
65 : : g_inputdeck.get< tag::ncomp >() ),
66 [ + - ]: 299 : m_npoin( Disc()->Coord()[0].size() ),
67 : : m_diag(),
68 : : m_stage( 0 ),
69 : : m_uc(),
70 : : m_pc(),
71 : : m_initial( 1 ),
72 : : m_uElemfields( m_u.nunk(), m_rhs.nprop() ),
73 : : m_pElemfields(m_u.nunk(),
74 : 299 : m_p.nprop()/g_inputdeck.get< tag::rdof >()),
75 : : m_uNodefields( m_npoin, m_rhs.nprop() ),
76 : : m_pNodefields(m_npoin,
77 : 299 : m_p.nprop()/g_inputdeck.get< tag::rdof >()),
78 : : m_uNodefieldsc(),
79 : : m_pNodefieldsc(),
80 : : m_boxelems(),
81 [ + - ]: 299 : m_srcFlag(m_u.nunk(), 0),
82 : : m_nrk(0),
83 : 598 : m_dte(m_u.nunk(), 0.0),
84 [ + - ][ + - ]: 1794 : m_finished(0)
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
85 : : // *****************************************************************************
86 : : // Constructor
87 : : //! \param[in] disc Discretization proxy
88 : : //! \param[in] bface Boundary-faces mapped to side set ids
89 : : //! \param[in] triinpoel Boundary-face connectivity
90 : : // *****************************************************************************
91 : : {
92 : : //! Runge-Kutta coefficients
93 [ + - ]: 299 : m_nrk = 2;
94 [ + - ]: 299 : m_rkcoef[0].resize(m_nrk);
95 [ + - ]: 299 : m_rkcoef[1].resize(m_nrk);
96 [ + - ]: 299 : if (m_nrk == 2) {
97 [ + - ][ + - ]: 598 : m_rkcoef = {{ {{ 0.0, 1.0/2.0 }}, {{ 1.0, 1.0/2.0 }} }};
98 : : }
99 : : else {
100 [ - - ][ - - ]: 0 : m_rkcoef = {{ {{ 0.0, 3.0/4.0, 1.0/3.0 }}, {{ 1.0, 1.0/4.0, 2.0/3.0 }} }};
101 : : }
102 : :
103 [ + - ]: 299 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
104 [ + + ]: 299 : g_inputdeck.get< tag::cmd, tag::quiescence >())
105 [ + - ][ + - ]: 306 : stateProxy.ckLocalBranch()->insert( "FV", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
[ - + ]
106 : : "FV" );
107 : :
108 : 299 : usesAtSync = true; // enable migration at AtSync
109 : :
110 : : // Enable SDAG wait for initially building the solution vector and limiting
111 [ + - ]: 299 : if (m_initial) {
112 [ + - ][ + - ]: 299 : thisProxy[ thisIndex ].wait4sol();
113 [ + - ][ + - ]: 299 : thisProxy[ thisIndex ].wait4lim();
114 [ + - ][ + - ]: 598 : thisProxy[ thisIndex ].wait4nod();
115 : : }
116 : :
117 [ + - ][ + - ]: 598 : m_ghosts[thisIndex].insert(m_disc, bface, triinpoel, m_u.nunk(),
118 [ + - ][ + - ]: 897 : CkCallback(CkIndex_FV::resizeSolVectors(), thisProxy[thisIndex]));
[ - + ][ - + ]
[ - - ][ - - ]
119 : :
120 : : // global-sync to call doneInserting on m_ghosts
121 [ + - ]: 299 : auto meshid = Disc()->MeshId();
122 [ + - ]: 299 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
123 [ + - ][ - - ]: 299 : CkCallback(CkReductionTarget(Transporter,doneInsertingGhosts),
124 [ + - ]: 299 : Disc()->Tr()) );
125 : 299 : }
126 : :
127 : : void
128 : 523 : FV::registerReducers()
129 : : // *****************************************************************************
130 : : // Configure Charm++ reduction types
131 : : //! \details Since this is a [initnode] routine, the runtime system executes the
132 : : //! routine exactly once on every logical node early on in the Charm++ init
133 : : //! sequence. Must be static as it is called without an object. See also:
134 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
135 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
136 : : // *****************************************************************************
137 : : {
138 : 523 : ElemDiagnostics::registerReducers();
139 : 523 : }
140 : :
141 : : void
142 : 385 : FV::ResumeFromSync()
143 : : // *****************************************************************************
144 : : // Return from migration
145 : : //! \details This is called when load balancing (LB) completes. The presence of
146 : : //! this function does not affect whether or not we block on LB.
147 : : // *****************************************************************************
148 : : {
149 [ - + ][ - - ]: 385 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
150 : :
151 [ + - ]: 385 : if (!g_inputdeck.get< tag::cmd, tag::nonblocking >()) next();
152 : 385 : }
153 : :
154 : : void
155 : 299 : FV::resizeSolVectors()
156 : : // *****************************************************************************
157 : : // Resize solution vectors after extension due to Ghosts
158 : : // *****************************************************************************
159 : : {
160 : : // Resize solution vectors, lhs and rhs by the number of ghost tets
161 : 299 : m_u.resize( myGhosts()->m_nunk );
162 : 299 : m_un.resize( myGhosts()->m_nunk );
163 : 299 : m_srcFlag.resize( myGhosts()->m_nunk );
164 : 299 : m_p.resize( myGhosts()->m_nunk );
165 : 299 : m_rhs.resize( myGhosts()->m_nunk );
166 : 299 : m_dte.resize( myGhosts()->m_nunk );
167 : :
168 : : // Size communication buffer for solution
169 [ + + ]: 897 : for (auto& u : m_uc) u.resize( myGhosts()->m_bid.size() );
170 [ + + ]: 897 : for (auto& p : m_pc) p.resize( myGhosts()->m_bid.size() );
171 : :
172 : : // Ensure that we also have all the geometry and connectivity data
173 : : // (including those of ghosts)
174 : : Assert( myGhosts()->m_geoElem.nunk() == m_u.nunk(),
175 : : "GeoElem unknowns size mismatch" );
176 : :
177 : : // Signal the runtime system that all workers have received their adjacency
178 : 299 : std::vector< std::size_t > meshdata{ myGhosts()->m_initial, Disc()->MeshId() };
179 : 299 : contribute( meshdata, CkReduction::sum_ulong,
180 [ + - ][ + - ]: 897 : CkCallback(CkReductionTarget(Transporter,comfinal), Disc()->Tr()) );
[ + - ][ - - ]
181 : 299 : }
182 : :
183 : : void
184 : 299 : FV::setup()
185 : : // *****************************************************************************
186 : : // Set initial conditions, generate lhs, output mesh
187 : : // *****************************************************************************
188 : : {
189 [ + - ]: 299 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
190 [ + + ]: 299 : g_inputdeck.get< tag::cmd, tag::quiescence >())
191 [ + - ][ + - ]: 306 : stateProxy.ckLocalBranch()->insert( "FV", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
[ - + ]
192 : : "setup" );
193 : :
194 : 299 : auto d = Disc();
195 : :
196 : : // Determine elements inside user-defined IC box
197 : 299 : g_fvpde[d->MeshId()].IcBoxElems( myGhosts()->m_geoElem,
198 : 299 : myGhosts()->m_fd.Esuel().size()/4, m_boxelems );
199 : :
200 : : // Compute volume of user-defined box IC
201 [ + - ][ - + ]: 299 : d->boxvol( {}, {}, 0 ); // punt for now
202 : :
203 : : // Query time history field output labels from all PDEs integrated
204 : : const auto& hist_points = g_inputdeck.get< tag::history_output, tag::point >();
205 [ - + ]: 299 : if (!hist_points.empty()) {
206 : 0 : std::vector< std::string > histnames;
207 [ - - ]: 0 : auto n = g_fvpde[d->MeshId()].histNames();
208 [ - - ]: 0 : histnames.insert( end(histnames), begin(n), end(n) );
209 [ - - ]: 0 : d->histheader( std::move(histnames) );
210 : : }
211 : 299 : }
212 : :
213 : : void
214 : 299 : FV::box( tk::real v, const std::vector< tk::real >& )
215 : : // *****************************************************************************
216 : : // Receive total box IC volume and set conditions in box
217 : : //! \param[in] v Total volume within user-specified box
218 : : // *****************************************************************************
219 : : {
220 : 299 : auto d = Disc();
221 : :
222 : : // Store user-defined box IC volume
223 : 299 : d->Boxvol() = v;
224 : :
225 : : // Set initial conditions for all PDEs
226 : 598 : g_fvpde[d->MeshId()].initialize( myGhosts()->m_geoElem, myGhosts()->m_inpoel,
227 : 299 : myGhosts()->m_coord, m_boxelems, d->ElemBlockId(), m_u, d->T(),
228 : 299 : myGhosts()->m_fd.Esuel().size()/4 );
229 : 299 : g_fvpde[d->MeshId()].updatePrimitives( m_u, m_p,
230 : 299 : myGhosts()->m_fd.Esuel().size()/4 );
231 : :
232 : : m_un = m_u;
233 : :
234 : : // Output initial conditions to file (regardless of whether it was requested)
235 [ + - ][ + - ]: 897 : startFieldOutput( CkCallback(CkIndex_FV::start(), thisProxy[thisIndex]) );
[ - + ][ - - ]
236 : 299 : }
237 : :
238 : : void
239 : 299 : FV::start()
240 : : // *****************************************************************************
241 : : // Start time stepping
242 : : // *****************************************************************************
243 : : {
244 : : // Start timer measuring time stepping wall clock time
245 : 299 : Disc()->Timer().zero();
246 : : // Zero grind-timer
247 : 299 : Disc()->grindZero();
248 : : // Start time stepping by computing the size of the next time step)
249 : 299 : next();
250 : 299 : }
251 : :
252 : : void
253 : 983 : FV::startFieldOutput( CkCallback c )
254 : : // *****************************************************************************
255 : : // Start preparing fields for output to file
256 : : //! \param[in] c Function to continue with after the write
257 : : // *****************************************************************************
258 : : {
259 : : // No field output in benchmark mode or if field output frequency not hit
260 [ + + ][ + + ]: 983 : if (g_inputdeck.get< tag::cmd, tag::benchmark >() || !fieldOutput()) {
261 : :
262 : 924 : c.send();
263 : :
264 : : } else {
265 : :
266 : : // Optionally refine mesh for field output
267 : 59 : auto d = Disc();
268 : :
269 [ - + ]: 59 : if (refinedOutput()) {
270 : :
271 : 0 : const auto& tr = tk::remap( myGhosts()->m_fd.Triinpoel(), d->Gid() );
272 [ - - ][ - - ]: 0 : d->Ref()->outref( myGhosts()->m_fd.Bface(), {}, tr, c );
[ - - ][ - - ]
[ - - ]
273 : :
274 : : } else {
275 : :
276 : : // cut off ghosts from mesh connectivity and coordinates
277 [ + - ][ - + ]: 118 : extractFieldOutput( {}, d->Chunk(), d->Coord(), {}, {}, d->NodeCommMap(),
[ - - ]
278 : : {}, {}, {}, c );
279 : :
280 : : }
281 : :
282 : : }
283 : 983 : }
284 : :
285 : : void
286 : 1368 : FV::next()
287 : : // *****************************************************************************
288 : : // Advance equations to next time step
289 : : // *****************************************************************************
290 : : {
291 : : // communicate solution ghost data (if any)
292 [ + + ]: 1368 : if (myGhosts()->m_sendGhost.empty())
293 : 128 : comsol_complete();
294 : : else
295 [ + + ]: 16592 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
296 [ + - ]: 14112 : std::vector< std::size_t > tetid( ghostdata.size() );
297 [ + - ][ + - ]: 28224 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
[ - - ]
298 [ + - ]: 14112 : prim( ghostdata.size() );
299 : : std::size_t j = 0;
300 [ + + ]: 174936 : for(const auto& i : ghostdata) {
301 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4,
302 : : "Sending solution ghost data" );
303 [ + - ]: 160824 : tetid[j] = i;
304 [ + - ][ - + ]: 160824 : u[j] = m_u[i];
305 [ + - ][ - + ]: 160824 : prim[j] = m_p[i];
306 : 160824 : ++j;
307 : : }
308 [ + - ][ + - ]: 28224 : thisProxy[ cid ].comsol( thisIndex, tetid, u, prim );
309 : : }
310 : :
311 : 1368 : ownsol_complete();
312 : 1368 : }
313 : :
314 : : void
315 : 14112 : FV::comsol( int fromch,
316 : : const std::vector< std::size_t >& tetid,
317 : : const std::vector< std::vector< tk::real > >& u,
318 : : const std::vector< std::vector< tk::real > >& prim )
319 : : // *****************************************************************************
320 : : // Receive chare-boundary solution ghost data from neighboring chares
321 : : //! \param[in] fromch Sender chare id
322 : : //! \param[in] tetid Ghost tet ids we receive solution data for
323 : : //! \param[in] u Solution ghost data
324 : : //! \param[in] prim Primitive variables in ghost cells
325 : : //! \details This function receives contributions to the unlimited solution
326 : : //! from fellow chares.
327 : : // *****************************************************************************
328 : : {
329 : : Assert( u.size() == tetid.size(), "Size mismatch in FV::comsol()" );
330 : : Assert( prim.size() == tetid.size(), "Size mismatch in FV::comsol()" );
331 : :
332 : : // Find local-to-ghost tet id map for sender chare
333 : 14112 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
334 : :
335 [ + + ]: 174936 : for (std::size_t i=0; i<tetid.size(); ++i) {
336 : 160824 : auto j = tk::cref_find( n, tetid[i] );
337 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
338 : : "Receiving solution non-ghost data" );
339 : 160824 : auto b = tk::cref_find( myGhosts()->m_bid, j );
340 : : Assert( b < m_uc[0].size(), "Indexing out of bounds" );
341 : 160824 : m_uc[0][b] = u[i];
342 : 160824 : m_pc[0][b] = prim[i];
343 : : }
344 : :
345 : : // if we have received all solution ghost contributions from neighboring
346 : : // chares (chares we communicate along chare-boundary faces with), and
347 : : // contributed our solution to these neighbors, proceed to reconstructions
348 [ + + ]: 14112 : if (++m_nsol == myGhosts()->m_sendGhost.size()) {
349 : 1240 : m_nsol = 0;
350 : 1240 : comsol_complete();
351 : : }
352 : 14112 : }
353 : :
354 : : void
355 [ + - ]: 59 : FV::extractFieldOutput(
356 : : const std::vector< std::size_t >& /*ginpoel*/,
357 : : const tk::UnsMesh::Chunk& chunk,
358 : : const tk::UnsMesh::Coords& coord,
359 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /*addedNodes*/,
360 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
361 : : const tk::NodeCommMap& nodeCommMap,
362 : : const std::map< int, std::vector< std::size_t > >& /* bface */,
363 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
364 : : const std::vector< std::size_t >& /* triinpoel */,
365 : : CkCallback c )
366 : : // *****************************************************************************
367 : : // Extract field output going to file
368 : : //! \param[in] chunk Field-output mesh chunk (connectivity and global<->local
369 : : //! id maps)
370 : : //! \param[in] coord Field-output mesh node coordinates
371 : : //! \param[in] addedTets Field-output mesh cells and their parents (local ids)
372 : : //! \param[in] nodeCommMap Field-output mesh node communication map
373 : : //! \param[in] c Function to continue with after the write
374 : : // *****************************************************************************
375 : : {
376 : : const auto& inpoel = std::get< 0 >( chunk );
377 : :
378 : : // Evaluate element solution on incoming mesh
379 [ + - ][ + - ]: 118 : evalSolution( *Disc(), inpoel, coord, addedTets, std::vector< std::size_t>{},
[ + + ]
380 [ + - ]: 59 : m_u, m_p, m_uElemfields, m_pElemfields, m_uNodefields, m_pNodefields );
381 : :
382 : : // Send node fields contributions to neighbor chares
383 [ + + ]: 59 : if (nodeCommMap.empty())
384 : 15 : comnodeout_complete();
385 : : else {
386 : : const auto& lid = std::get< 2 >( chunk );
387 : 88 : auto esup = tk::genEsup( inpoel, 4 );
388 [ + + ]: 176 : for(const auto& [ch,nodes] : nodeCommMap) {
389 : : // Pack node field data in chare boundary nodes
390 : : std::vector< std::vector< tk::real > >
391 [ + - ][ + - ]: 396 : lu( m_uNodefields.nprop(), std::vector< tk::real >( nodes.size() ) );
[ + - ]
392 : : std::vector< std::vector< tk::real > >
393 [ + - ][ + - ]: 396 : lp( m_pNodefields.nprop(), std::vector< tk::real >( nodes.size() ) );
394 [ + + ]: 1320 : for (std::size_t f=0; f<m_uNodefields.nprop(); ++f) {
395 : : std::size_t j = 0;
396 [ + + ]: 26532 : for (auto g : nodes)
397 : 25344 : lu[f][j++] = m_uNodefields(tk::cref_find(lid,g),f);
398 : : }
399 [ + + ]: 792 : for (std::size_t f=0; f<m_pNodefields.nprop(); ++f) {
400 : : std::size_t j = 0;
401 [ + + ]: 14740 : for (auto g : nodes)
402 : 14080 : lp[f][j++] = m_pNodefields(tk::cref_find(lid,g),f);
403 : : }
404 : : // Pack (partial) number of elements surrounding chare boundary nodes
405 [ + - ]: 132 : std::vector< std::size_t > nesup( nodes.size() );
406 : : std::size_t j = 0;
407 [ + + ]: 2948 : for (auto g : nodes) {
408 : 2816 : auto i = tk::cref_find( lid, g );
409 : 2816 : nesup[j++] = esup.second[i+1] - esup.second[i];
410 : : }
411 [ + - ][ + - ]: 396 : thisProxy[ch].comnodeout(
[ + - ][ - - ]
412 [ + - ][ - + ]: 264 : std::vector<std::size_t>(begin(nodes),end(nodes)), nesup, lu, lp );
413 : : }
414 : : }
415 : :
416 [ + - ]: 59 : ownnod_complete( c );
417 : 59 : }
418 : :
419 : : void
420 : 1368 : FV::reco()
421 : : // *****************************************************************************
422 : : // Compute reconstructions
423 : : // *****************************************************************************
424 : : {
425 : 1368 : const auto rdof = g_inputdeck.get< tag::rdof >();
426 : :
427 : : // Combine own and communicated contributions of unreconstructed solution and
428 : : // degrees of freedom in cells (if p-adaptive)
429 [ + + ]: 163560 : for (const auto& b : myGhosts()->m_bid) {
430 : : Assert( m_uc[0][b.second].size() == m_u.nprop(), "ncomp size mismatch" );
431 : : Assert( m_pc[0][b.second].size() == m_p.nprop(), "ncomp size mismatch" );
432 [ + + ]: 7661496 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
433 : 7500672 : m_u(b.first,c) = m_uc[0][b.second][c];
434 : : }
435 [ + + ]: 3947640 : for (std::size_t c=0; c<m_p.nprop(); ++c) {
436 : 3786816 : m_p(b.first,c) = m_pc[0][b.second][c];
437 : : }
438 : : }
439 : :
440 [ + - ]: 1368 : if (rdof > 1) {
441 : : // Reconstruct second-order solution and primitive quantities
442 : 2736 : g_fvpde[Disc()->MeshId()].reconstruct( myGhosts()->m_geoElem, myGhosts()->m_fd,
443 : 1368 : myGhosts()->m_esup, myGhosts()->m_inpoel, myGhosts()->m_coord, m_u, m_p );
444 : : }
445 : :
446 : : // start limiting
447 : 1368 : lim();
448 : 1368 : }
449 : :
450 : : void
451 : 1368 : FV::lim()
452 : : // *****************************************************************************
453 : : // Compute limiter function
454 : : // *****************************************************************************
455 : : {
456 : 1368 : const auto rdof = g_inputdeck.get< tag::rdof >();
457 : :
458 [ + - ]: 1368 : if (rdof > 1) {
459 : 2736 : g_fvpde[Disc()->MeshId()].limit( myGhosts()->m_geoFace, myGhosts()->m_fd,
460 : 1368 : myGhosts()->m_esup, myGhosts()->m_inpoel, m_srcFlag, m_u, m_p );
461 : : }
462 : :
463 : : // Send limited solution to neighboring chares
464 [ + + ]: 1368 : if (myGhosts()->m_sendGhost.empty())
465 : 128 : comlim_complete();
466 : : else
467 [ + + ]: 16592 : for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
468 [ + - ]: 14112 : std::vector< std::size_t > tetid( ghostdata.size() );
469 [ + - ][ + - ]: 28224 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
[ - - ]
470 [ + - ]: 14112 : prim( ghostdata.size() );
471 : : std::size_t j = 0;
472 [ + + ]: 174936 : for(const auto& i : ghostdata) {
473 : : Assert( i < myGhosts()->m_fd.Esuel().size()/4,
474 : : "Sending limiter ghost data" );
475 [ + - ]: 160824 : tetid[j] = i;
476 [ + - ][ - + ]: 160824 : u[j] = m_u[i];
477 [ + - ][ - + ]: 160824 : prim[j] = m_p[i];
478 : 160824 : ++j;
479 : : }
480 [ + - ][ + - ]: 28224 : thisProxy[ cid ].comlim( thisIndex, tetid, u, prim );
481 : : }
482 : :
483 : 1368 : ownlim_complete();
484 : 1368 : }
485 : :
486 : : void
487 : 14112 : FV::comlim( int fromch,
488 : : const std::vector< std::size_t >& tetid,
489 : : const std::vector< std::vector< tk::real > >& u,
490 : : const std::vector< std::vector< tk::real > >& prim )
491 : : // *****************************************************************************
492 : : // Receive chare-boundary limiter ghost data from neighboring chares
493 : : //! \param[in] fromch Sender chare id
494 : : //! \param[in] tetid Ghost tet ids we receive solution data for
495 : : //! \param[in] u Limited high-order solution
496 : : //! \param[in] prim Limited high-order primitive quantities
497 : : //! \details This function receives contributions to the limited solution from
498 : : //! fellow chares.
499 : : // *****************************************************************************
500 : : {
501 : : Assert( u.size() == tetid.size(), "Size mismatch in FV::comlim()" );
502 : : Assert( prim.size() == tetid.size(), "Size mismatch in FV::comlim()" );
503 : :
504 : : // Find local-to-ghost tet id map for sender chare
505 : 14112 : const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
506 : :
507 [ + + ]: 174936 : for (std::size_t i=0; i<tetid.size(); ++i) {
508 : 160824 : auto j = tk::cref_find( n, tetid[i] );
509 : : Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
510 : : "Receiving solution non-ghost data" );
511 : 160824 : auto b = tk::cref_find( myGhosts()->m_bid, j );
512 : : Assert( b < m_uc[1].size(), "Indexing out of bounds" );
513 : : Assert( b < m_pc[1].size(), "Indexing out of bounds" );
514 : 160824 : m_uc[1][b] = u[i];
515 : 160824 : m_pc[1][b] = prim[i];
516 : : }
517 : :
518 : : // if we have received all solution ghost contributions from neighboring
519 : : // chares (chares we communicate along chare-boundary faces with), and
520 : : // contributed our solution to these neighbors, proceed to limiting
521 [ + + ]: 14112 : if (++m_nlim == myGhosts()->m_sendGhost.size()) {
522 : 1240 : m_nlim = 0;
523 : 1240 : comlim_complete();
524 : : }
525 : 14112 : }
526 : :
527 : : void
528 : 1368 : FV::dt()
529 : : // *****************************************************************************
530 : : // Compute time step size
531 : : // *****************************************************************************
532 : : {
533 : 1368 : auto d = Disc();
534 : :
535 : : // Combine own and communicated contributions of limited solution and degrees
536 : : // of freedom in cells (if p-adaptive)
537 [ + + ]: 163560 : for (const auto& b : myGhosts()->m_bid) {
538 : : Assert( m_uc[1][b.second].size() == m_u.nprop(), "ncomp size mismatch" );
539 : : Assert( m_pc[1][b.second].size() == m_p.nprop(), "ncomp size mismatch" );
540 [ + + ]: 7661496 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
541 : 7500672 : m_u(b.first,c) = m_uc[1][b.second][c];
542 : : }
543 [ + + ]: 3947640 : for (std::size_t c=0; c<m_p.nprop(); ++c) {
544 : 3786816 : m_p(b.first,c) = m_pc[1][b.second][c];
545 : : }
546 : : }
547 : :
548 : 1368 : auto mindt = std::numeric_limits< tk::real >::max();
549 : :
550 [ + + ]: 1368 : if (m_stage == 0)
551 : : {
552 [ + + ]: 684 : auto const_dt = g_inputdeck.get< tag::dt >();
553 : : auto eps = std::numeric_limits< tk::real >::epsilon();
554 : :
555 : : // use constant dt if configured
556 [ + + ]: 684 : if (std::abs(const_dt) > eps) {
557 : :
558 : 584 : mindt = const_dt;
559 : :
560 : : } else { // compute dt based on CFL
561 : :
562 : : // find the minimum dt across all PDEs integrated
563 : : auto eqdt =
564 : 200 : g_fvpde[d->MeshId()].dt( myGhosts()->m_fd, myGhosts()->m_geoFace,
565 : 100 : myGhosts()->m_geoElem, m_u, m_p, myGhosts()->m_fd.Esuel().size()/4,
566 : 100 : m_srcFlag, m_dte );
567 [ + - ]: 100 : if (eqdt < mindt) mindt = eqdt;
568 : :
569 : : // time-step suppression for unsteady problems
570 : : tk::real coeff(1.0);
571 [ + - ]: 100 : if (!g_inputdeck.get< tag::steady_state >()) {
572 : 100 : auto ramp_steps = g_inputdeck.get< tag::cfl_ramping_steps >();
573 [ + - ][ + - ]: 100 : if (g_inputdeck.get< tag::cfl_ramping >() && d->It() < ramp_steps)
574 : 100 : coeff = 1.0/static_cast< tk::real >(ramp_steps)
575 : 100 : * static_cast< tk::real >(d->It());
576 : : }
577 : : else {
578 [ - - ]: 0 : for (auto& edt : m_dte) edt *= g_inputdeck.get< tag::cfl >();
579 : : }
580 : :
581 : 100 : mindt *= coeff * g_inputdeck.get< tag::cfl >();
582 : : }
583 : : }
584 : : else
585 : : {
586 : 684 : mindt = d->Dt();
587 : : }
588 : :
589 : : // Contribute to minimum dt across all chares then advance to next step
590 [ + - ]: 1368 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
591 : 1368 : CkCallback(CkReductionTarget(FV,solve), thisProxy) );
592 : 1368 : }
593 : :
594 : : void
595 : 1368 : FV::solve( tk::real newdt )
596 : : // *****************************************************************************
597 : : // Compute right-hand side of discrete transport equations
598 : : //! \param[in] newdt Size of this new time step
599 : : // *****************************************************************************
600 : : {
601 : : // Enable SDAG wait for building the solution vector during the next stage
602 [ + - ]: 1368 : thisProxy[ thisIndex ].wait4sol();
603 [ + - ]: 1368 : thisProxy[ thisIndex ].wait4lim();
604 [ + - ]: 1368 : thisProxy[ thisIndex ].wait4nod();
605 : :
606 : 1368 : auto d = Disc();
607 : 1368 : const auto rdof = g_inputdeck.get< tag::rdof >();
608 : 1368 : const auto neq = m_u.nprop()/rdof;
609 : :
610 : : // Set new time step size
611 [ + + ]: 1368 : if (m_stage == 0) d->setdt( newdt );
612 : :
613 : : // Update Un
614 [ + + ]: 1368 : if (m_stage == 0) m_un = m_u;
615 : :
616 : : // physical time at time-stage for computing exact source terms for
617 : : // unsteady problems
618 : 1368 : tk::real physT(d->T());
619 : : // 2-stage RK
620 [ + - ]: 1368 : if (m_nrk == 2) {
621 [ + + ]: 1368 : if (m_stage == 1) {
622 : 684 : physT += d->Dt();
623 : : }
624 : : }
625 : : // 3-stage RK
626 : : else {
627 [ - - ]: 0 : if (m_stage == 1) {
628 : 0 : physT += d->Dt();
629 : : }
630 [ - - ]: 0 : else if (m_stage == 2) {
631 : 0 : physT += 0.5*d->Dt();
632 : : }
633 : : }
634 : :
635 : : // Compute rhs
636 : 2736 : g_fvpde[d->MeshId()].rhs( physT, myGhosts()->m_geoFace, myGhosts()->m_geoElem,
637 : 1368 : myGhosts()->m_fd, myGhosts()->m_inpoel, myGhosts()->m_coord,
638 : 1368 : d->ElemBlockId(), m_u, m_p, m_rhs, m_srcFlag );
639 : :
640 : : // Explicit time-stepping using RK3 to discretize time-derivative
641 : 1368 : const auto steady = g_inputdeck.get< tag::steady_state >();
642 [ + + ]: 426992 : for (std::size_t e=0; e<myGhosts()->m_nunk; ++e) {
643 : 425624 : auto vole = myGhosts()->m_geoElem(e,0);
644 [ + + ]: 4841672 : for (std::size_t c=0; c<neq; ++c)
645 : : {
646 : 4416048 : auto dte = d->Dt();
647 [ - + ]: 4416048 : if (steady) dte = m_dte[e];
648 : 4416048 : auto rmark = c*rdof;
649 [ + - ]: 4416048 : m_u(e, rmark) = m_rkcoef[0][m_stage] * m_un(e, rmark)
650 : 4416048 : + m_rkcoef[1][m_stage] * ( m_u(e, rmark)
651 : 4416048 : + dte * m_rhs(e, c)/vole );
652 : : // zero out reconstructed dofs of equations using reduced dofs
653 [ + - ]: 4416048 : if (rdof > 1) {
654 [ + + ]: 17664192 : for (std::size_t k=1; k<rdof; ++k)
655 : : {
656 : 13248144 : rmark = c*rdof+k;
657 : 13248144 : m_u(e, rmark) = 0.0;
658 : : }
659 : : }
660 : : }
661 : : }
662 : :
663 : : // Update primitives based on the evolved solution
664 : 1368 : g_fvpde[d->MeshId()].updatePrimitives( m_u, m_p,
665 : 1368 : myGhosts()->m_fd.Esuel().size()/4 );
666 [ + - ]: 1368 : if (!g_inputdeck.get< tag::accuracy_test >()) {
667 : 1368 : g_fvpde[d->MeshId()].cleanTraceMaterial( physT, myGhosts()->m_geoElem, m_u,
668 : 1368 : m_p, myGhosts()->m_fd.Esuel().size()/4 );
669 : : }
670 : :
671 [ + + ]: 1368 : if (m_stage < m_nrk-1) {
672 : :
673 : : // continue with next time step stage
674 : 684 : stage();
675 : :
676 : : } else {
677 : :
678 : : // Increase number of iterations and physical time
679 : 684 : d->next();
680 : :
681 : : // Compute diagnostics, e.g., residuals
682 [ + - ]: 684 : auto diag_computed = m_diag.compute( *d,
683 [ + - ][ + - ]: 684 : m_u.nunk()-myGhosts()->m_fd.Esuel().size()/4, myGhosts()->m_geoElem,
[ + - ]
684 [ + - ]: 684 : std::vector< std::size_t>{}, m_u, m_un );
685 : :
686 : : // Continue to mesh refinement (if configured)
687 [ - + ][ - - ]: 684 : if (!diag_computed) refine( std::vector< tk::real >( m_u.nprop(), 1.0 ) );
[ - - ]
688 : :
689 : : }
690 : 1368 : }
691 : :
692 : : void
693 : 684 : FV::refine( const std::vector< tk::real >& l2res )
694 : : // *****************************************************************************
695 : : // Optionally refine/derefine mesh
696 : : //! \param[in] l2res L2-norms of the residual for each scalar component
697 : : //! computed across the whole problem
698 : : // *****************************************************************************
699 : : {
700 : 684 : auto d = Disc();
701 : :
702 : : // Assess convergence for steady state
703 : 684 : const auto steady = g_inputdeck.get< tag::steady_state >();
704 : 684 : const auto residual = g_inputdeck.get< tag::residual >();
705 : 684 : const auto rc = g_inputdeck.get< tag::rescomp >() - 1;
706 : :
707 : : bool converged(false);
708 [ - + ]: 684 : if (steady) converged = l2res[rc] < residual;
709 : :
710 : : // this is the last time step if max time of max number of time steps
711 : : // reached or the residual has reached its convergence criterion
712 [ + + ][ - + ]: 684 : if (d->finished() or converged) m_finished = 1;
713 : :
714 : 684 : auto dtref = g_inputdeck.get< tag::amr, tag::dtref >();
715 : 684 : auto dtfreq = g_inputdeck.get< tag::amr, tag::dtfreq >();
716 : :
717 : : // if t>0 refinement enabled and we hit the dtref frequency
718 [ - + ][ - - ]: 684 : if (dtref && !(d->It() % dtfreq)) { // refine
719 : :
720 : 0 : d->startvol();
721 [ - - ][ - - ]: 0 : d->Ref()->dtref( myGhosts()->m_fd.Bface(), {},
[ - - ][ - - ]
722 : 0 : tk::remap(myGhosts()->m_fd.Triinpoel(),d->Gid()) );
723 : 0 : d->refined() = 1;
724 : :
725 : : } else { // do not refine
726 : :
727 : 684 : d->refined() = 0;
728 : 684 : stage();
729 : :
730 : : }
731 : 684 : }
732 : :
733 : : void
734 : 0 : FV::resizePostAMR(
735 : : const std::vector< std::size_t >& /*ginpoel*/,
736 : : const tk::UnsMesh::Chunk& chunk,
737 : : const tk::UnsMesh::Coords& coord,
738 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /*addedNodes*/,
739 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
740 : : const std::set< std::size_t >& removedNodes,
741 : : const std::unordered_map< std::size_t, std::size_t >& amrNodeMap,
742 : : const tk::NodeCommMap& nodeCommMap,
743 : : const std::map< int, std::vector< std::size_t > >& bface,
744 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
745 : : const std::vector< std::size_t >& triinpoel,
746 : : const std::unordered_map< std::size_t, std::set< std::size_t > >& elemblockid )
747 : : // *****************************************************************************
748 : : // Receive new mesh from Refiner
749 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
750 : : //! \param[in] coord New mesh node coordinates
751 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
752 : : //! \param[in] removedNodes Newly removed mesh node local ids
753 : : //! \param[in] amrNodeMap Node id map after amr (local ids)
754 : : //! \param[in] nodeCommMap New node communication map
755 : : //! \param[in] bface Boundary-faces mapped to side set ids
756 : : //! \param[in] triinpoel Boundary-face connectivity
757 : : //! \param[in] elemblockid Local tet ids associated with mesh block ids
758 : : // *****************************************************************************
759 : : {
760 : 0 : auto d = Disc();
761 : :
762 : : // Set flag that indicates that we are during time stepping
763 : 0 : m_initial = 0;
764 : 0 : myGhosts()->m_initial = 0;
765 : :
766 : : // Zero field output iteration count between two mesh refinement steps
767 : 0 : d->Itf() = 0;
768 : :
769 : : // Increase number of iterations with mesh refinement
770 : 0 : ++d->Itr();
771 : :
772 : : // Save old number of elements
773 : 0 : [[maybe_unused]] auto old_nelem = myGhosts()->m_inpoel.size()/4;
774 : :
775 : : // Resize mesh data structures
776 : 0 : d->resizePostAMR( chunk, coord, amrNodeMap, nodeCommMap, removedNodes,
777 : : elemblockid );
778 : :
779 : : // Update state
780 : 0 : myGhosts()->m_inpoel = d->Inpoel();
781 : 0 : myGhosts()->m_coord = d->Coord();
782 : 0 : auto nelem = myGhosts()->m_inpoel.size()/4;
783 : : m_p.resize( nelem );
784 : : m_u.resize( nelem );
785 : 0 : m_srcFlag.resize( nelem );
786 : : m_un.resize( nelem );
787 : : m_rhs.resize( nelem );
788 : :
789 [ - - ][ - - ]: 0 : myGhosts()->m_fd = FaceData( myGhosts()->m_inpoel, bface,
[ - - ][ - - ]
[ - - ]
790 : 0 : tk::remap(triinpoel,d->Lid()) );
791 : :
792 [ - - ]: 0 : myGhosts()->m_geoFace =
793 : 0 : tk::Fields( tk::genGeoFaceTri( myGhosts()->m_fd.Nipfac(),
794 : 0 : myGhosts()->m_fd.Inpofa(), coord ) );
795 [ - - ]: 0 : myGhosts()->m_geoElem = tk::Fields( tk::genGeoElemTet( myGhosts()->m_inpoel,
796 : 0 : coord ) );
797 : :
798 : 0 : myGhosts()->m_nfac = myGhosts()->m_fd.Inpofa().size()/3;
799 : 0 : myGhosts()->m_nunk = nelem;
800 : 0 : m_npoin = coord[0].size();
801 : 0 : myGhosts()->m_bndFace.clear();
802 : 0 : myGhosts()->m_exptGhost.clear();
803 : 0 : myGhosts()->m_sendGhost.clear();
804 : 0 : myGhosts()->m_ghost.clear();
805 : 0 : myGhosts()->m_esup.clear();
806 : :
807 : : // Update solution on new mesh, P0 (cell center value) only for now
808 : : m_un = m_u;
809 : : auto pn = m_p;
810 : 0 : auto unprop = m_u.nprop();
811 : : auto pnprop = m_p.nprop();
812 [ - - ]: 0 : for (const auto& [child,parent] : addedTets) {
813 : : Assert( child < nelem, "Indexing out of new solution vector" );
814 : : Assert( parent < old_nelem, "Indexing out of old solution vector" );
815 [ - - ]: 0 : for (std::size_t i=0; i<unprop; ++i) m_u(child,i) = m_un(parent,i);
816 [ - - ]: 0 : for (std::size_t i=0; i<pnprop; ++i) m_p(child,i) = pn(parent,i);
817 : : }
818 : : m_un = m_u;
819 : :
820 : : // Resize communication buffers
821 [ - - ][ - - ]: 0 : m_ghosts[thisIndex].resizeComm();
[ - - ][ - - ]
822 : 0 : }
823 : :
824 : : bool
825 : 107 : FV::fieldOutput() const
826 : : // *****************************************************************************
827 : : // Decide wether to output field data
828 : : //! \return True if field data is output in this step
829 : : // *****************************************************************************
830 : : {
831 : 107 : auto d = Disc();
832 : :
833 : : // Output field data
834 [ + + ][ + - ]: 107 : return d->fielditer() or d->fieldtime() or d->fieldrange() or m_finished;
[ + - ][ - + ]
835 : : }
836 : :
837 : : bool
838 : 59 : FV::refinedOutput() const
839 : : // *****************************************************************************
840 : : // Decide if we write field output using a refined mesh
841 : : //! \return True if field output will use a refined mesh
842 : : // *****************************************************************************
843 : : {
844 [ - + ]: 59 : return g_inputdeck.get< tag::field_output, tag::refined >() &&
845 [ - - ]: 0 : g_inputdeck.get< tag::scheme >() != ctr::SchemeType::FV;
846 : : }
847 : :
848 : : void
849 : 59 : FV::writeFields( CkCallback c )
850 : : // *****************************************************************************
851 : : // Output mesh field data
852 : : //! \param[in] c Function to continue with after the write
853 : : // *****************************************************************************
854 : : {
855 : 59 : auto d = Disc();
856 : :
857 : : const auto& inpoel = std::get< 0 >( d->Chunk() );
858 : 118 : auto esup = tk::genEsup( inpoel, 4 );
859 : 59 : auto nelem = inpoel.size() / 4;
860 : :
861 : : // Combine own and communicated contributions and finish averaging of node
862 : : // field output in chare boundary nodes
863 : : const auto& lid = std::get< 2 >( d->Chunk() );
864 [ + + ]: 2611 : for (const auto& [g,f] : m_uNodefieldsc) {
865 : : Assert( m_uNodefields.nprop() == f.first.size(), "Size mismatch" );
866 : 2552 : auto p = tk::cref_find( lid, g );
867 [ + + ]: 25520 : for (std::size_t i=0; i<f.first.size(); ++i) {
868 : 22968 : m_uNodefields(p,i) += f.first[i];
869 : 22968 : m_uNodefields(p,i) /= static_cast< tk::real >(
870 : 22968 : esup.second[p+1] - esup.second[p] + f.second );
871 : : }
872 : : }
873 : 59 : tk::destroy( m_uNodefieldsc );
874 [ + + ]: 2611 : for (const auto& [g,f] : m_pNodefieldsc) {
875 : : Assert( m_pNodefields.nprop() == f.first.size(), "Size mismatch" );
876 : 2552 : auto p = tk::cref_find( lid, g );
877 [ + + ]: 15312 : for (std::size_t i=0; i<f.first.size(); ++i) {
878 : 12760 : m_pNodefields(p,i) += f.first[i];
879 : 12760 : m_pNodefields(p,i) /= static_cast< tk::real >(
880 : 12760 : esup.second[p+1] - esup.second[p] + f.second );
881 : : }
882 : : }
883 : 59 : tk::destroy( m_pNodefieldsc );
884 : :
885 : : // Lambda to decide if a node (global id) is on a chare boundary of the field
886 : : // output mesh. p - global node id, return true if node is on the chare
887 : : // boundary.
888 : 14075 : auto chbnd = [ this ]( std::size_t p ) {
889 : : return
890 : 14075 : std::any_of( Disc()->NodeCommMap().cbegin(), Disc()->NodeCommMap().cend(),
891 : 14075 : [&](const auto& s) { return s.second.find(p) != s.second.cend(); } );
892 : 59 : };
893 : :
894 : : // Finish computing node field output averages in internal nodes
895 : 59 : auto npoin = d->Coord()[0].size();
896 : : auto& gid = std::get< 1 >( d->Chunk() );
897 [ + + ]: 14134 : for (std::size_t p=0; p<npoin; ++p) {
898 [ + - ][ + + ]: 14075 : if (!chbnd(gid[p])) {
899 : 11523 : auto n = static_cast< tk::real >( esup.second[p+1] - esup.second[p] );
900 [ + + ]: 115230 : for (std::size_t i=0; i<m_uNodefields.nprop(); ++i)
901 : 103707 : m_uNodefields(p,i) /= n;
902 [ + + ]: 69138 : for (std::size_t i=0; i<m_pNodefields.nprop(); ++i)
903 : 57615 : m_pNodefields(p,i) /= n;
904 : : }
905 : : }
906 : :
907 : : // Collect field output from numerical solution requested by user
908 : 59 : auto elemfields = numericFieldOutput( m_uElemfields, tk::Centering::ELEM,
909 [ + - ][ + - ]: 118 : g_fvpde[Disc()->MeshId()].OutVarFn(), m_pElemfields );
[ + - ]
910 : 59 : auto nodefields = numericFieldOutput( m_uNodefields, tk::Centering::NODE,
911 [ + - ][ + - ]: 177 : g_fvpde[Disc()->MeshId()].OutVarFn(), m_pNodefields );
[ + - ][ + - ]
912 : :
913 : : // Collect field output from analytical solutions (if exist)
914 : : const auto& coord = d->Coord();
915 [ + - ]: 59 : auto geoElem = tk::genGeoElemTet( inpoel, coord );
916 [ + - ]: 59 : auto t = Disc()->T();
917 [ + - ]: 59 : analyticFieldOutput( g_fvpde[d->MeshId()], tk::Centering::ELEM,
918 [ + - ][ + - ]: 236 : geoElem.extract_comp(1), geoElem.extract_comp(2), geoElem.extract_comp(3),
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
919 : : t, elemfields );
920 [ + - ]: 59 : analyticFieldOutput( g_fvpde[d->MeshId()], tk::Centering::NODE, coord[0],
921 : : coord[1], coord[2], t, nodefields );
922 : :
923 : : // Add sound speed vector
924 [ + - ][ - - ]: 59 : std::vector< tk::real > soundspd(nelem, 0.0);
925 [ + - ]: 59 : g_fvpde[d->MeshId()].soundspeed(nelem, m_u, m_p, soundspd);
926 [ + - ]: 59 : elemfields.push_back(soundspd);
927 : :
928 : : // Add source flag array to element-centered field output
929 [ + - ][ - - ]: 59 : std::vector< tk::real > srcFlag( begin(m_srcFlag), end(m_srcFlag) );
930 : : // Here m_srcFlag has a size of m_u.nunk() which is the number of the
931 : : // elements within this partition (nelem) plus the ghost partition cells.
932 : : // For the purpose of output, we only need the solution data within this
933 : : // partition. Therefore, resizing it to nelem removes the extra partition
934 : : // boundary allocations in the srcFlag vector. Since the code assumes that
935 : : // the boundary elements are on the top, the resize operation keeps the lower
936 : : // portion.
937 [ + - ]: 59 : srcFlag.resize( nelem );
938 [ + - ]: 59 : elemfields.push_back( srcFlag );
939 : :
940 : : // Query fields names requested by user
941 [ + - ]: 118 : auto elemfieldnames = numericFieldNames( tk::Centering::ELEM );
942 [ + - ]: 118 : auto nodefieldnames = numericFieldNames( tk::Centering::NODE );
943 : :
944 : : // Collect field output names for analytical solutions
945 [ + - ]: 59 : analyticFieldNames( g_fvpde[d->MeshId()], tk::Centering::ELEM, elemfieldnames );
946 [ + - ]: 59 : analyticFieldNames( g_fvpde[d->MeshId()], tk::Centering::NODE, nodefieldnames );
947 : :
948 [ + - ]: 59 : elemfieldnames.push_back( "sound speed" );
949 [ + - ]: 59 : elemfieldnames.push_back( "src_flag" );
950 : :
951 : : Assert( elemfieldnames.size() == elemfields.size(), "Size mismatch" );
952 : : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
953 : :
954 : : // Collect surface output names
955 [ + - ]: 118 : auto surfnames = g_fvpde[d->MeshId()].surfNames();
956 : :
957 : : // Collect surface field solution
958 [ + - ]: 59 : const auto& fd = myGhosts()->m_fd;
959 [ + - ]: 118 : auto elemsurfs = g_fvpde[d->MeshId()].surfOutput(fd, m_u, m_p);
960 : :
961 : : // Output chare mesh and fields metadata to file
962 [ + - ]: 59 : const auto& triinpoel = tk::remap( fd.Triinpoel(), d->Gid() );
963 [ + - ][ + - ]: 177 : d->write( inpoel, d->Coord(), fd.Bface(), {},
[ + - ][ + - ]
[ - - ][ - - ]
964 [ + - ]: 118 : tk::remap( triinpoel, lid ), elemfieldnames, nodefieldnames,
965 : : surfnames, {}, elemfields, nodefields, elemsurfs, {}, c );
966 : 59 : }
967 : :
968 : : void
969 : 132 : FV::comnodeout( const std::vector< std::size_t >& gid,
970 : : const std::vector< std::size_t >& nesup,
971 : : const std::vector< std::vector< tk::real > >& Lu,
972 : : const std::vector< std::vector< tk::real > >& Lp )
973 : : // *****************************************************************************
974 : : // Receive chare-boundary nodal solution (for field output) contributions from
975 : : // neighboring chares
976 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
977 : : //! \param[in] nesup Number of elements surrounding points
978 : : //! \param[in] Lu Partial contributions of solution nodal fields to
979 : : //! chare-boundary nodes
980 : : //! \param[in] Lp Partial contributions of primitive quantity nodal fields to
981 : : //! chare-boundary nodes
982 : : // *****************************************************************************
983 : : {
984 : : Assert( gid.size() == nesup.size(), "Size mismatch" );
985 : : Assert(Lu.size() == m_uNodefields.nprop(), "Fields size mismatch");
986 : : Assert(Lp.size() == m_pNodefields.nprop(), "Fields size mismatch");
987 [ + + ]: 1320 : for (std::size_t f=0; f<Lu.size(); ++f)
988 : : Assert( gid.size() == Lu[f].size(), "Size mismatch" );
989 [ + + ]: 792 : for (std::size_t f=0; f<Lp.size(); ++f)
990 : : Assert( gid.size() == Lp[f].size(), "Size mismatch" );
991 : :
992 [ + + ]: 2948 : for (std::size_t i=0; i<gid.size(); ++i) {
993 : : auto& nfu = m_uNodefieldsc[ gid[i] ];
994 : 2816 : nfu.first.resize( Lu.size() );
995 [ + + ]: 28160 : for (std::size_t f=0; f<Lu.size(); ++f) nfu.first[f] += Lu[f][i];
996 : 2816 : nfu.second += nesup[i];
997 : 2816 : auto& nfp = m_pNodefieldsc[ gid[i] ];
998 : 2816 : nfp.first.resize( Lp.size() );
999 [ + + ]: 16896 : for (std::size_t f=0; f<Lp.size(); ++f) nfp.first[f] += Lp[f][i];
1000 : 2816 : nfp.second += nesup[i];
1001 : : }
1002 : :
1003 : : // When we have heard from all chares we communicate with, this chare is done
1004 [ + + ]: 132 : if (++m_nnod == Disc()->NodeCommMap().size()) {
1005 : 44 : m_nnod = 0;
1006 : 44 : comnodeout_complete();
1007 : : }
1008 : 132 : }
1009 : :
1010 : : void
1011 : 1368 : FV::stage()
1012 : : // *****************************************************************************
1013 : : // Evaluate whether to continue with next time step stage
1014 : : // *****************************************************************************
1015 : : {
1016 : : // Increment Runge-Kutta stage counter
1017 : 1368 : ++m_stage;
1018 : :
1019 : : // if not all Runge-Kutta stages complete, continue to next time stage,
1020 : : // otherwise prepare for nodal field output
1021 [ + + ]: 1368 : if (m_stage < m_nrk)
1022 : 684 : next();
1023 : : else
1024 [ + - ][ + - ]: 2052 : startFieldOutput( CkCallback(CkIndex_FV::step(), thisProxy[thisIndex]) );
[ - + ][ - - ]
1025 : 1368 : }
1026 : :
1027 : : void
1028 : 385 : FV::evalLB( int nrestart )
1029 : : // *****************************************************************************
1030 : : // Evaluate whether to do load balancing
1031 : : //! \param[in] nrestart Number of times restarted
1032 : : // *****************************************************************************
1033 : : {
1034 : 385 : auto d = Disc();
1035 : :
1036 : : // Detect if just returned from a checkpoint and if so, zero timers and flag
1037 [ + + ]: 385 : if (d->restarted( nrestart )) m_finished = 0;
1038 : :
1039 : 385 : const auto lbfreq = g_inputdeck.get< tag::cmd, tag::lbfreq >();
1040 : 385 : const auto nonblocking = g_inputdeck.get< tag::cmd, tag::nonblocking >();
1041 : :
1042 : : // Load balancing if user frequency is reached or after the second time-step
1043 [ - + ][ - - ]: 385 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
1044 : :
1045 : 385 : AtSync();
1046 [ - + ]: 385 : if (nonblocking) next();
1047 : :
1048 : : } else {
1049 : :
1050 : 0 : next();
1051 : :
1052 : : }
1053 : 385 : }
1054 : :
1055 : : void
1056 : 380 : FV::evalRestart()
1057 : : // *****************************************************************************
1058 : : // Evaluate whether to save checkpoint/restart
1059 : : // *****************************************************************************
1060 : : {
1061 : 380 : auto d = Disc();
1062 : :
1063 : 380 : const auto rsfreq = g_inputdeck.get< tag::cmd, tag::rsfreq >();
1064 : 380 : const auto benchmark = g_inputdeck.get< tag::cmd, tag::benchmark >();
1065 : :
1066 [ + + ][ - + ]: 380 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
1067 : :
1068 : 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
1069 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
1070 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
[ - - ]
1071 : :
1072 : : } else {
1073 : :
1074 : 380 : evalLB( /* nrestart = */ -1 );
1075 : :
1076 : : }
1077 : 380 : }
1078 : :
1079 : : void
1080 : 684 : FV::step()
1081 : : // *****************************************************************************
1082 : : // Evaluate wether to continue with next time step
1083 : : // *****************************************************************************
1084 : : {
1085 : 684 : auto d = Disc();
1086 : :
1087 : : // Output time history
1088 [ + - ][ + - ]: 684 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ]
1089 : 0 : std::vector< std::vector< tk::real > > hist;
1090 [ - - ][ - - ]: 0 : auto h = g_fvpde[d->MeshId()].histOutput( d->Hist(), myGhosts()->m_inpoel,
1091 [ - - ][ - - ]: 0 : myGhosts()->m_coord, m_u, m_p );
1092 [ - - ]: 0 : hist.insert( end(hist), begin(h), end(h) );
1093 [ - - ]: 0 : d->history( std::move(hist) );
1094 : : }
1095 : :
1096 : : // Output one-liner status report to screen
1097 : 684 : d->status();
1098 : : // Reset Runge-Kutta stage counter
1099 : 684 : m_stage = 0;
1100 : :
1101 : : // If neither max iterations nor max time reached, continue, otherwise finish
1102 [ + + ]: 684 : if (not m_finished) {
1103 : :
1104 : 380 : evalRestart();
1105 : :
1106 : : } else {
1107 : :
1108 : 304 : auto meshid = d->MeshId();
1109 [ + - ]: 608 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1110 : 608 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
1111 : :
1112 : : }
1113 : 684 : }
1114 : :
1115 : : #include "NoWarning/fv.def.h"
|
