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