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 : 298 : 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 : 298 : 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 [ + - ]: 298 : m_u( Disc()->Inpoel().size()/4,
59 : 298 : g_inputdeck.get< tag::rdof >()*
60 : 298 : g_inputdeck.get< tag::ncomp >() ),
61 : : m_un( m_u.nunk(), m_u.nprop() ),
62 : 298 : m_p( m_u.nunk(), g_inputdeck.get< tag::rdof >()*
63 [ + - ][ + - ]: 298 : 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 [ + - ]: 298 : 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 : 298 : m_p.nprop()/g_inputdeck.get< tag::rdof >()),
76 : : m_uNodefields( m_npoin, m_lhs.nprop() ),
77 : : m_pNodefields(m_npoin,
78 : 298 : m_p.nprop()/g_inputdeck.get< tag::rdof >()),
79 : : m_uNodefieldsc(),
80 : : m_pNodefieldsc(),
81 : : m_boxelems(),
82 [ + - ]: 298 : m_srcFlag(m_u.nunk(), 0),
83 : : m_nrk(0),
84 : 596 : m_dte(m_u.nunk(), 0.0),
85 [ + - ][ + - ]: 1788 : 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 [ + - ]: 298 : m_nrk = 2;
95 [ + - ]: 298 : m_rkcoef[0].resize(m_nrk);
96 [ + - ]: 298 : m_rkcoef[1].resize(m_nrk);
97 [ + - ]: 298 : if (m_nrk == 2) {
98 [ + - ][ + - ]: 596 : 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 [ + - ]: 298 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
105 [ + + ]: 298 : g_inputdeck.get< tag::cmd, tag::quiescence >())
106 [ + - ][ + - ]: 304 : stateProxy.ckLocalBranch()->insert( "FV", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
[ - + ]
107 : : "FV" );
108 : :
109 : 298 : usesAtSync = true; // enable migration at AtSync
110 : :
111 : : // Enable SDAG wait for initially building the solution vector and limiting
112 [ + - ]: 298 : if (m_initial) {
113 [ + - ][ + - ]: 298 : thisProxy[ thisIndex ].wait4sol();
114 [ + - ][ + - ]: 298 : thisProxy[ thisIndex ].wait4lim();
115 [ + - ][ + - ]: 596 : thisProxy[ thisIndex ].wait4nod();
116 : : }
117 : :
118 [ + - ][ + - ]: 596 : m_ghosts[thisIndex].insert(m_disc, bface, triinpoel, m_u.nunk(),
119 [ + - ][ + - ]: 894 : CkCallback(CkIndex_FV::resizeSolVectors(), thisProxy[thisIndex]));
[ - + ][ - + ]
[ - - ][ - - ]
120 : :
121 : : // global-sync to call doneInserting on m_ghosts
122 [ + - ]: 298 : auto meshid = Disc()->MeshId();
123 [ + - ]: 298 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
124 [ + - ][ - - ]: 298 : CkCallback(CkReductionTarget(Transporter,doneInsertingGhosts),
125 [ + - ]: 298 : Disc()->Tr()) );
126 : 298 : }
127 : :
128 : : void
129 : 572 : 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 : 572 : ElemDiagnostics::registerReducers();
140 : 572 : }
141 : :
142 : : void
143 : 536 : 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 [ - + ][ - - ]: 536 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
151 : :
152 [ + - ]: 536 : if (!g_inputdeck.get< tag::cmd, tag::nonblocking >()) next();
153 : 536 : }
154 : :
155 : : void
156 : 298 : 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 : 298 : m_u.resize( myGhosts()->m_nunk );
163 : 298 : m_un.resize( myGhosts()->m_nunk );
164 : 298 : m_srcFlag.resize( myGhosts()->m_nunk );
165 : 298 : m_p.resize( myGhosts()->m_nunk );
166 : 298 : m_lhs.resize( myGhosts()->m_nunk );
167 : 298 : m_rhs.resize( myGhosts()->m_nunk );
168 : 298 : m_dte.resize( myGhosts()->m_nunk );
169 : :
170 : : // Size communication buffer for solution
171 [ + + ]: 894 : for (auto& u : m_uc) u.resize( myGhosts()->m_bid.size() );
172 [ + + ]: 894 : 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 : 298 : std::vector< std::size_t > meshdata{ myGhosts()->m_initial, Disc()->MeshId() };
181 : 298 : contribute( meshdata, CkReduction::sum_ulong,
182 [ + - ][ + - ]: 894 : CkCallback(CkReductionTarget(Transporter,comfinal), Disc()->Tr()) );
[ + - ][ - - ]
183 : 298 : }
184 : :
185 : : void
186 : 298 : FV::setup()
187 : : // *****************************************************************************
188 : : // Set initial conditions, generate lhs, output mesh
189 : : // *****************************************************************************
190 : : {
191 [ + - ]: 298 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
192 [ + + ]: 298 : g_inputdeck.get< tag::cmd, tag::quiescence >())
193 [ + - ][ + - ]: 304 : stateProxy.ckLocalBranch()->insert( "FV", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
[ - + ]
194 : : "setup" );
195 : :
196 : 298 : 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 : 298 : lhs();
204 : :
205 : : // Determine elements inside user-defined IC box
206 : 298 : g_fvpde[d->MeshId()].IcBoxElems( myGhosts()->m_geoElem,
207 : 298 : myGhosts()->m_fd.Esuel().size()/4, m_boxelems );
208 : :
209 : : // Compute volume of user-defined box IC
210 [ + - ][ - + ]: 298 : 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 [ - + ]: 298 : 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 : 298 : }
221 : :
222 : : void
223 : 298 : 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 : 298 : auto d = Disc();
230 : :
231 : : // Store user-defined box IC volume
232 : 298 : d->Boxvol() = v;
233 : :
234 : : // Set initial conditions for all PDEs
235 : 596 : g_fvpde[d->MeshId()].initialize( m_lhs, myGhosts()->m_inpoel,
236 : 298 : myGhosts()->m_coord, m_boxelems, d->ElemBlockId(), m_u, d->T(),
237 : 298 : myGhosts()->m_fd.Esuel().size()/4 );
238 : 298 : g_fvpde[d->MeshId()].updatePrimitives( m_u, m_p,
239 : 298 : 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 [ + - ][ + - ]: 894 : startFieldOutput( CkCallback(CkIndex_FV::start(), thisProxy[thisIndex]) );
[ - + ][ - - ]
245 : 298 : }
246 : :
247 : : void
248 : 298 : FV::start()
249 : : // *****************************************************************************
250 : : // Start time stepping
251 : : // *****************************************************************************
252 : : {
253 : : // Start timer measuring time stepping wall clock time
254 : 298 : Disc()->Timer().zero();
255 : : // Zero grind-timer
256 : 298 : Disc()->grindZero();
257 : : // Start time stepping by computing the size of the next time step)
258 : 298 : next();
259 : 298 : }
260 : :
261 : : void
262 : 1132 : 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 [ + + ][ + + ]: 1132 : if (g_inputdeck.get< tag::cmd, tag::benchmark >() || !fieldOutput()) {
270 : :
271 : 1104 : c.send();
272 : :
273 : : } else {
274 : :
275 : : // Optionally refine mesh for field output
276 : 28 : auto d = Disc();
277 : :
278 [ - + ]: 28 : 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 [ + - ][ - + ]: 56 : extractFieldOutput( {}, d->Chunk(), d->Coord(), {}, {}, d->NodeCommMap(),
[ - - ]
287 : : {}, {}, {}, c );
288 : :
289 : : }
290 : :
291 : : }
292 : 1132 : }
293 : :
294 : : void
295 : 1668 : FV::next()
296 : : // *****************************************************************************
297 : : // Advance equations to next time step
298 : : // *****************************************************************************
299 : : {
300 : : // communicate solution ghost data (if any)
301 [ + + ]: 1668 : if (myGhosts()->m_sendGhost.empty())
302 : 108 : 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 : 1668 : ownsol_complete();
321 : 1668 : }
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 [ + - ]: 28 : 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 [ + - ][ + - ]: 56 : evalSolution( *Disc(), inpoel, coord, addedTets, std::vector< std::size_t>{},
[ + + ]
389 [ + - ]: 28 : m_u, m_p, m_uElemfields, m_pElemfields, m_uNodefields, m_pNodefields );
390 : :
391 : : // Send node fields contributions to neighbor chares
392 [ + + ]: 28 : if (nodeCommMap.empty())
393 : 4 : 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 [ + - ]: 28 : ownnod_complete( c );
426 : 28 : }
427 : :
428 : : void
429 : 298 : FV::lhs()
430 : : // *****************************************************************************
431 : : // Compute left-hand side of discrete transport equations
432 : : // *****************************************************************************
433 : : {
434 : 298 : g_fvpde[Disc()->MeshId()].lhs( myGhosts()->m_geoElem, m_lhs );
435 : :
436 [ - + ]: 298 : if (!m_initial) stage();
437 : 298 : }
438 : :
439 : : void
440 : 1668 : FV::reco()
441 : : // *****************************************************************************
442 : : // Compute reconstructions
443 : : // *****************************************************************************
444 : : {
445 : 1668 : 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 [ + + ]: 351220 : 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 [ + - ]: 1668 : if (rdof > 1) {
461 : : // Reconstruct second-order solution and primitive quantities
462 : 3336 : g_fvpde[Disc()->MeshId()].reconstruct( myGhosts()->m_geoElem, myGhosts()->m_fd,
463 : 1668 : myGhosts()->m_esup, myGhosts()->m_inpoel, myGhosts()->m_coord, m_u, m_p );
464 : : }
465 : :
466 : : // start limiting
467 : 1668 : lim();
468 : 1668 : }
469 : :
470 : : void
471 : 1668 : FV::lim()
472 : : // *****************************************************************************
473 : : // Compute limiter function
474 : : // *****************************************************************************
475 : : {
476 : 1668 : const auto rdof = g_inputdeck.get< tag::rdof >();
477 : :
478 [ + - ]: 1668 : if (rdof > 1) {
479 : 3336 : g_fvpde[Disc()->MeshId()].limit( myGhosts()->m_geoFace, myGhosts()->m_fd,
480 : 1668 : myGhosts()->m_esup,
481 : 1668 : myGhosts()->m_inpoel, myGhosts()->m_coord, m_srcFlag, m_u, m_p );
482 : : }
483 : :
484 : : // Send limited solution to neighboring chares
485 [ + + ]: 1668 : if (myGhosts()->m_sendGhost.empty())
486 : 108 : 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 : 1668 : ownlim_complete();
505 : 1668 : }
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 : 1668 : FV::dt()
550 : : // *****************************************************************************
551 : : // Compute time step size
552 : : // *****************************************************************************
553 : : {
554 : 1668 : auto d = Disc();
555 : :
556 : : // Combine own and communicated contributions of limited solution and degrees
557 : : // of freedom in cells (if p-adaptive)
558 [ + + ]: 351220 : 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 : 1668 : auto mindt = std::numeric_limits< tk::real >::max();
570 : :
571 [ + + ]: 1668 : if (m_stage == 0)
572 : : {
573 [ + + ]: 834 : 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 [ + + ]: 834 : 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 : 500 : g_fvpde[d->MeshId()].dt( myGhosts()->m_fd, myGhosts()->m_geoFace,
586 : 250 : myGhosts()->m_geoElem, m_u, m_p, myGhosts()->m_fd.Esuel().size()/4,
587 : 250 : m_srcFlag, m_dte );
588 [ + - ]: 250 : if (eqdt < mindt) mindt = eqdt;
589 : :
590 : : // time-step suppression for unsteady problems
591 : : tk::real coeff(1.0);
592 [ + - ]: 250 : if (!g_inputdeck.get< tag::steady_state >()) {
593 [ + - ]: 250 : if (d->It() < 100) coeff = 0.01 * static_cast< tk::real >(d->It());
594 : : }
595 : : else {
596 [ - - ]: 0 : for (auto& edt : m_dte) edt *= g_inputdeck.get< tag::cfl >();
597 : : }
598 : :
599 : 250 : mindt *= coeff * g_inputdeck.get< tag::cfl >();
600 : : }
601 : : }
602 : : else
603 : : {
604 : 834 : mindt = d->Dt();
605 : : }
606 : :
607 : : // Contribute to minimum dt across all chares then advance to next step
608 [ + - ]: 1668 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
609 : 1668 : CkCallback(CkReductionTarget(FV,solve), thisProxy) );
610 : 1668 : }
611 : :
612 : : void
613 : 1668 : FV::solve( tk::real newdt )
614 : : // *****************************************************************************
615 : : // Compute right-hand side of discrete transport equations
616 : : //! \param[in] newdt Size of this new time step
617 : : // *****************************************************************************
618 : : {
619 : : // Enable SDAG wait for building the solution vector during the next stage
620 [ + - ]: 1668 : thisProxy[ thisIndex ].wait4sol();
621 [ + - ]: 1668 : thisProxy[ thisIndex ].wait4lim();
622 [ + - ]: 1668 : thisProxy[ thisIndex ].wait4nod();
623 : :
624 : 1668 : auto d = Disc();
625 : 1668 : const auto rdof = g_inputdeck.get< tag::rdof >();
626 : 1668 : const auto neq = m_u.nprop()/rdof;
627 : :
628 : : // Set new time step size
629 [ + + ]: 1668 : if (m_stage == 0) d->setdt( newdt );
630 : :
631 : : // Update Un
632 [ + + ]: 1668 : if (m_stage == 0) m_un = m_u;
633 : :
634 : : // physical time at time-stage for computing exact source terms for
635 : : // unsteady problems
636 : 1668 : tk::real physT(d->T());
637 : : // 2-stage RK
638 [ + - ]: 1668 : if (m_nrk == 2) {
639 [ + + ]: 1668 : if (m_stage == 1) {
640 : 834 : physT += d->Dt();
641 : : }
642 : : }
643 : : // 3-stage RK
644 : : else {
645 [ - - ]: 0 : if (m_stage == 1) {
646 : 0 : physT += d->Dt();
647 : : }
648 [ - - ]: 0 : else if (m_stage == 2) {
649 : 0 : physT += 0.5*d->Dt();
650 : : }
651 : : }
652 : :
653 : : // Compute rhs
654 : 3336 : g_fvpde[d->MeshId()].rhs( physT, myGhosts()->m_geoFace, myGhosts()->m_geoElem,
655 : 1668 : myGhosts()->m_fd, myGhosts()->m_inpoel, myGhosts()->m_coord,
656 : 1668 : d->ElemBlockId(), m_u, m_p, m_rhs, m_srcFlag );
657 : :
658 : : // Explicit time-stepping using RK3 to discretize time-derivative
659 : 1668 : const auto steady = g_inputdeck.get< tag::steady_state >();
660 [ + + ]: 646112 : for (std::size_t e=0; e<myGhosts()->m_nunk; ++e)
661 [ + + ]: 7029872 : for (std::size_t c=0; c<neq; ++c)
662 : : {
663 : 6385428 : auto dte = d->Dt();
664 [ - + ]: 6385428 : if (steady) dte = m_dte[e];
665 : 6385428 : auto rmark = c*rdof;
666 [ + - ]: 6385428 : m_u(e, rmark) = m_rkcoef[0][m_stage] * m_un(e, rmark)
667 : 6385428 : + m_rkcoef[1][m_stage] * ( m_u(e, rmark)
668 : 6385428 : + dte * m_rhs(e, c)/m_lhs(e, c) );
669 : : // zero out reconstructed dofs of equations using reduced dofs
670 [ + - ]: 6385428 : if (rdof > 1) {
671 [ + + ]: 25541712 : for (std::size_t k=1; k<rdof; ++k)
672 : : {
673 : 19156284 : rmark = c*rdof+k;
674 : 19156284 : m_u(e, rmark) = 0.0;
675 : : }
676 : : }
677 : : }
678 : :
679 : : // Update primitives based on the evolved solution
680 : 1668 : g_fvpde[d->MeshId()].updatePrimitives( m_u, m_p,
681 : 1668 : myGhosts()->m_fd.Esuel().size()/4 );
682 [ + - ]: 1668 : if (!g_inputdeck.get< tag::accuracy_test >()) {
683 : 1668 : g_fvpde[d->MeshId()].cleanTraceMaterial( physT, myGhosts()->m_geoElem, m_u,
684 : 1668 : m_p, myGhosts()->m_fd.Esuel().size()/4 );
685 : : }
686 : :
687 [ + + ]: 1668 : if (m_stage < m_nrk-1) {
688 : :
689 : : // continue with next time step stage
690 : 834 : stage();
691 : :
692 : : } else {
693 : :
694 : : // Increase number of iterations and physical time
695 : 834 : d->next();
696 : :
697 : : // Compute diagnostics, e.g., residuals
698 [ + - ]: 834 : auto diag_computed = m_diag.compute( *d,
699 [ + - ][ + - ]: 834 : m_u.nunk()-myGhosts()->m_fd.Esuel().size()/4, myGhosts()->m_geoElem,
[ + - ]
700 [ + - ]: 834 : std::vector< std::size_t>{}, m_u, m_un );
701 : :
702 : : // Continue to mesh refinement (if configured)
703 [ + + ][ + - ]: 990 : if (!diag_computed) refine( std::vector< tk::real >( m_u.nprop(), 1.0 ) );
[ + - ]
704 : :
705 : : }
706 : 1668 : }
707 : :
708 : : void
709 : 834 : FV::refine( const std::vector< tk::real >& l2res )
710 : : // *****************************************************************************
711 : : // Optionally refine/derefine mesh
712 : : //! \param[in] l2res L2-norms of the residual for each scalar component
713 : : //! computed across the whole problem
714 : : // *****************************************************************************
715 : : {
716 : 834 : auto d = Disc();
717 : :
718 : : // Assess convergence for steady state
719 : 834 : const auto steady = g_inputdeck.get< tag::steady_state >();
720 : 834 : const auto residual = g_inputdeck.get< tag::residual >();
721 : 834 : const auto rc = g_inputdeck.get< tag::rescomp >() - 1;
722 : :
723 : : bool converged(false);
724 [ - + ]: 834 : if (steady) converged = l2res[rc] < residual;
725 : :
726 : : // this is the last time step if max time of max number of time steps
727 : : // reached or the residual has reached its convergence criterion
728 [ + + ][ - + ]: 834 : if (d->finished() or converged) m_finished = 1;
729 : :
730 : 834 : auto dtref = g_inputdeck.get< tag::amr, tag::dtref >();
731 : 834 : auto dtfreq = g_inputdeck.get< tag::amr, tag::dtfreq >();
732 : :
733 : : // if t>0 refinement enabled and we hit the dtref frequency
734 [ - + ][ - - ]: 834 : if (dtref && !(d->It() % dtfreq)) { // refine
735 : :
736 : 0 : d->startvol();
737 [ - - ][ - - ]: 0 : d->Ref()->dtref( myGhosts()->m_fd.Bface(), {},
[ - - ][ - - ]
738 : 0 : tk::remap(myGhosts()->m_fd.Triinpoel(),d->Gid()) );
739 : 0 : d->refined() = 1;
740 : :
741 : : } else { // do not refine
742 : :
743 : 834 : d->refined() = 0;
744 : 834 : stage();
745 : :
746 : : }
747 : 834 : }
748 : :
749 : : void
750 : 0 : FV::resizePostAMR(
751 : : const std::vector< std::size_t >& /*ginpoel*/,
752 : : const tk::UnsMesh::Chunk& chunk,
753 : : const tk::UnsMesh::Coords& coord,
754 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /*addedNodes*/,
755 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
756 : : const std::set< std::size_t >& removedNodes,
757 : : const std::unordered_map< std::size_t, std::size_t >& amrNodeMap,
758 : : const tk::NodeCommMap& nodeCommMap,
759 : : const std::map< int, std::vector< std::size_t > >& bface,
760 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
761 : : const std::vector< std::size_t >& triinpoel,
762 : : const std::unordered_map< std::size_t, std::set< std::size_t > >& elemblockid )
763 : : // *****************************************************************************
764 : : // Receive new mesh from Refiner
765 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
766 : : //! \param[in] coord New mesh node coordinates
767 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
768 : : //! \param[in] removedNodes Newly removed mesh node local ids
769 : : //! \param[in] amrNodeMap Node id map after amr (local ids)
770 : : //! \param[in] nodeCommMap New node communication map
771 : : //! \param[in] bface Boundary-faces mapped to side set ids
772 : : //! \param[in] triinpoel Boundary-face connectivity
773 : : //! \param[in] elemblockid Local tet ids associated with mesh block ids
774 : : // *****************************************************************************
775 : : {
776 : 0 : auto d = Disc();
777 : :
778 : : // Set flag that indicates that we are during time stepping
779 : 0 : m_initial = 0;
780 : 0 : myGhosts()->m_initial = 0;
781 : :
782 : : // Zero field output iteration count between two mesh refinement steps
783 : 0 : d->Itf() = 0;
784 : :
785 : : // Increase number of iterations with mesh refinement
786 : 0 : ++d->Itr();
787 : :
788 : : // Save old number of elements
789 : 0 : [[maybe_unused]] auto old_nelem = myGhosts()->m_inpoel.size()/4;
790 : :
791 : : // Resize mesh data structures
792 : 0 : d->resizePostAMR( chunk, coord, amrNodeMap, nodeCommMap, removedNodes,
793 : : elemblockid );
794 : :
795 : : // Update state
796 : 0 : myGhosts()->m_inpoel = d->Inpoel();
797 : 0 : myGhosts()->m_coord = d->Coord();
798 : 0 : auto nelem = myGhosts()->m_inpoel.size()/4;
799 : : m_p.resize( nelem );
800 : : m_u.resize( nelem );
801 : 0 : m_srcFlag.resize( nelem );
802 : : m_un.resize( nelem );
803 : : m_lhs.resize( nelem );
804 : : m_rhs.resize( nelem );
805 : :
806 [ - - ][ - - ]: 0 : myGhosts()->m_fd = FaceData( myGhosts()->m_inpoel, bface,
[ - - ][ - - ]
[ - - ]
807 : 0 : tk::remap(triinpoel,d->Lid()) );
808 : :
809 [ - - ]: 0 : myGhosts()->m_geoFace =
810 : 0 : tk::Fields( tk::genGeoFaceTri( myGhosts()->m_fd.Nipfac(),
811 : 0 : myGhosts()->m_fd.Inpofa(), coord ) );
812 [ - - ]: 0 : myGhosts()->m_geoElem = tk::Fields( tk::genGeoElemTet( myGhosts()->m_inpoel,
813 : 0 : coord ) );
814 : :
815 : 0 : myGhosts()->m_nfac = myGhosts()->m_fd.Inpofa().size()/3;
816 : 0 : myGhosts()->m_nunk = nelem;
817 : 0 : m_npoin = coord[0].size();
818 : 0 : myGhosts()->m_bndFace.clear();
819 : 0 : myGhosts()->m_exptGhost.clear();
820 : 0 : myGhosts()->m_sendGhost.clear();
821 : 0 : myGhosts()->m_ghost.clear();
822 : 0 : myGhosts()->m_esup.clear();
823 : :
824 : : // Update solution on new mesh, P0 (cell center value) only for now
825 : : m_un = m_u;
826 : : auto pn = m_p;
827 : 0 : auto unprop = m_u.nprop();
828 : : auto pnprop = m_p.nprop();
829 [ - - ]: 0 : for (const auto& [child,parent] : addedTets) {
830 : : Assert( child < nelem, "Indexing out of new solution vector" );
831 : : Assert( parent < old_nelem, "Indexing out of old solution vector" );
832 [ - - ]: 0 : for (std::size_t i=0; i<unprop; ++i) m_u(child,i) = m_un(parent,i);
833 [ - - ]: 0 : for (std::size_t i=0; i<pnprop; ++i) m_p(child,i) = pn(parent,i);
834 : : }
835 : : m_un = m_u;
836 : :
837 : : // Resize communication buffers
838 [ - - ][ - - ]: 0 : m_ghosts[thisIndex].resizeComm();
[ - - ][ - - ]
839 : 0 : }
840 : :
841 : : bool
842 : 256 : FV::fieldOutput() const
843 : : // *****************************************************************************
844 : : // Decide wether to output field data
845 : : //! \return True if field data is output in this step
846 : : // *****************************************************************************
847 : : {
848 : 256 : auto d = Disc();
849 : :
850 : : // Output field data
851 [ + + ][ + - ]: 256 : return d->fielditer() or d->fieldtime() or d->fieldrange() or m_finished;
[ + - ][ - + ]
852 : : }
853 : :
854 : : bool
855 : 28 : FV::refinedOutput() const
856 : : // *****************************************************************************
857 : : // Decide if we write field output using a refined mesh
858 : : //! \return True if field output will use a refined mesh
859 : : // *****************************************************************************
860 : : {
861 [ - + ]: 28 : return g_inputdeck.get< tag::field_output, tag::refined >() &&
862 [ - - ]: 0 : g_inputdeck.get< tag::scheme >() != ctr::SchemeType::FV;
863 : : }
864 : :
865 : : void
866 : 28 : FV::writeFields( CkCallback c )
867 : : // *****************************************************************************
868 : : // Output mesh field data
869 : : //! \param[in] c Function to continue with after the write
870 : : // *****************************************************************************
871 : : {
872 : 28 : auto d = Disc();
873 : :
874 : : const auto& inpoel = std::get< 0 >( d->Chunk() );
875 : 56 : auto esup = tk::genEsup( inpoel, 4 );
876 : 28 : auto nelem = inpoel.size() / 4;
877 : :
878 : : // Combine own and communicated contributions and finish averaging of node
879 : : // field output in chare boundary nodes
880 : : const auto& lid = std::get< 2 >( d->Chunk() );
881 [ + + ]: 2740 : for (const auto& [g,f] : m_uNodefieldsc) {
882 : : Assert( m_uNodefields.nprop() == f.first.size(), "Size mismatch" );
883 : 2712 : auto p = tk::cref_find( lid, g );
884 [ + + ]: 27120 : for (std::size_t i=0; i<f.first.size(); ++i) {
885 : 24408 : m_uNodefields(p,i) += f.first[i];
886 : 24408 : m_uNodefields(p,i) /= static_cast< tk::real >(
887 : 24408 : esup.second[p+1] - esup.second[p] + f.second );
888 : : }
889 : : }
890 : 28 : tk::destroy( m_uNodefieldsc );
891 [ + + ]: 2740 : for (const auto& [g,f] : m_pNodefieldsc) {
892 : : Assert( m_pNodefields.nprop() == f.first.size(), "Size mismatch" );
893 : 2712 : auto p = tk::cref_find( lid, g );
894 [ + + ]: 16272 : for (std::size_t i=0; i<f.first.size(); ++i) {
895 : 13560 : m_pNodefields(p,i) += f.first[i];
896 : 13560 : m_pNodefields(p,i) /= static_cast< tk::real >(
897 : 13560 : esup.second[p+1] - esup.second[p] + f.second );
898 : : }
899 : : }
900 : 28 : tk::destroy( m_pNodefieldsc );
901 : :
902 : : // Lambda to decide if a node (global id) is on a chare boundary of the field
903 : : // output mesh. p - global node id, return true if node is on the chare
904 : : // boundary.
905 : 5462 : auto chbnd = [ this ]( std::size_t p ) {
906 : : return
907 : 5462 : std::any_of( Disc()->NodeCommMap().cbegin(), Disc()->NodeCommMap().cend(),
908 : 5462 : [&](const auto& s) { return s.second.find(p) != s.second.cend(); } );
909 : 28 : };
910 : :
911 : : // Finish computing node field output averages in internal nodes
912 : 28 : auto npoin = d->Coord()[0].size();
913 : : auto& gid = std::get< 1 >( d->Chunk() );
914 [ + + ]: 5490 : for (std::size_t p=0; p<npoin; ++p) {
915 [ + - ][ + + ]: 5462 : if (!chbnd(gid[p])) {
916 : 2750 : auto n = static_cast< tk::real >( esup.second[p+1] - esup.second[p] );
917 [ + + ]: 27500 : for (std::size_t i=0; i<m_uNodefields.nprop(); ++i)
918 : 24750 : m_uNodefields(p,i) /= n;
919 [ + + ]: 16500 : for (std::size_t i=0; i<m_pNodefields.nprop(); ++i)
920 : 13750 : m_pNodefields(p,i) /= n;
921 : : }
922 : : }
923 : :
924 : : // Collect field output from numerical solution requested by user
925 : 28 : auto elemfields = numericFieldOutput( m_uElemfields, tk::Centering::ELEM,
926 [ + - ][ + - ]: 56 : g_fvpde[Disc()->MeshId()].OutVarFn(), m_pElemfields );
[ + - ]
927 : 28 : auto nodefields = numericFieldOutput( m_uNodefields, tk::Centering::NODE,
928 [ + - ][ + - ]: 84 : g_fvpde[Disc()->MeshId()].OutVarFn(), m_pNodefields );
[ + - ][ + - ]
929 : :
930 : : // Collect field output from analytical solutions (if exist)
931 : : const auto& coord = d->Coord();
932 [ + - ]: 28 : auto geoElem = tk::genGeoElemTet( inpoel, coord );
933 [ + - ]: 28 : auto t = Disc()->T();
934 [ + - ]: 28 : analyticFieldOutput( g_fvpde[d->MeshId()], tk::Centering::ELEM,
935 [ + - ][ + - ]: 112 : geoElem.extract_comp(1), geoElem.extract_comp(2), geoElem.extract_comp(3),
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
936 : : t, elemfields );
937 [ + - ]: 28 : analyticFieldOutput( g_fvpde[d->MeshId()], tk::Centering::NODE, coord[0],
938 : : coord[1], coord[2], t, nodefields );
939 : :
940 : : // Add sound speed vector
941 [ + - ][ - - ]: 28 : std::vector< tk::real > soundspd(nelem, 0.0);
942 [ + - ]: 28 : g_fvpde[d->MeshId()].soundspeed(nelem, m_u, m_p, soundspd);
943 [ + - ]: 28 : elemfields.push_back(soundspd);
944 : :
945 : : // Add source flag array to element-centered field output
946 [ + - ][ - - ]: 28 : std::vector< tk::real > srcFlag( begin(m_srcFlag), end(m_srcFlag) );
947 : : // Here m_srcFlag has a size of m_u.nunk() which is the number of the
948 : : // elements within this partition (nelem) plus the ghost partition cells.
949 : : // For the purpose of output, we only need the solution data within this
950 : : // partition. Therefore, resizing it to nelem removes the extra partition
951 : : // boundary allocations in the srcFlag vector. Since the code assumes that
952 : : // the boundary elements are on the top, the resize operation keeps the lower
953 : : // portion.
954 [ + - ]: 28 : srcFlag.resize( nelem );
955 [ + - ]: 28 : elemfields.push_back( srcFlag );
956 : :
957 : : // Query fields names requested by user
958 [ + - ]: 56 : auto elemfieldnames = numericFieldNames( tk::Centering::ELEM );
959 [ + - ]: 56 : auto nodefieldnames = numericFieldNames( tk::Centering::NODE );
960 : :
961 : : // Collect field output names for analytical solutions
962 [ + - ]: 28 : analyticFieldNames( g_fvpde[d->MeshId()], tk::Centering::ELEM, elemfieldnames );
963 [ + - ]: 28 : analyticFieldNames( g_fvpde[d->MeshId()], tk::Centering::NODE, nodefieldnames );
964 : :
965 [ + - ]: 28 : elemfieldnames.push_back( "sound speed" );
966 [ + - ]: 28 : elemfieldnames.push_back( "src_flag" );
967 : :
968 : : Assert( elemfieldnames.size() == elemfields.size(), "Size mismatch" );
969 : : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
970 : :
971 : : // Collect surface output names
972 [ + - ]: 56 : auto surfnames = g_fvpde[d->MeshId()].surfNames();
973 : :
974 : : // Collect surface field solution
975 [ + - ]: 28 : const auto& fd = myGhosts()->m_fd;
976 [ + - ]: 56 : auto elemsurfs = g_fvpde[d->MeshId()].surfOutput(fd, m_u, m_p);
977 : :
978 : : // Output chare mesh and fields metadata to file
979 [ + - ]: 28 : const auto& triinpoel = tk::remap( fd.Triinpoel(), d->Gid() );
980 [ + - ][ + - ]: 84 : d->write( inpoel, d->Coord(), fd.Bface(), {},
[ + - ][ + - ]
[ - - ][ - - ]
981 [ + - ]: 56 : tk::remap( triinpoel, lid ), elemfieldnames, nodefieldnames,
982 : : surfnames, {}, elemfields, nodefields, elemsurfs, {}, c );
983 : 28 : }
984 : :
985 : : void
986 : 72 : FV::comnodeout( const std::vector< std::size_t >& gid,
987 : : const std::vector< std::size_t >& nesup,
988 : : const std::vector< std::vector< tk::real > >& Lu,
989 : : const std::vector< std::vector< tk::real > >& Lp )
990 : : // *****************************************************************************
991 : : // Receive chare-boundary nodal solution (for field output) contributions from
992 : : // neighboring chares
993 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
994 : : //! \param[in] nesup Number of elements surrounding points
995 : : //! \param[in] Lu Partial contributions of solution nodal fields to
996 : : //! chare-boundary nodes
997 : : //! \param[in] Lp Partial contributions of primitive quantity nodal fields to
998 : : //! chare-boundary nodes
999 : : // *****************************************************************************
1000 : : {
1001 : : Assert( gid.size() == nesup.size(), "Size mismatch" );
1002 : : Assert(Lu.size() == m_uNodefields.nprop(), "Fields size mismatch");
1003 : : Assert(Lp.size() == m_pNodefields.nprop(), "Fields size mismatch");
1004 [ + + ]: 720 : for (std::size_t f=0; f<Lu.size(); ++f)
1005 : : Assert( gid.size() == Lu[f].size(), "Size mismatch" );
1006 [ + + ]: 432 : for (std::size_t f=0; f<Lp.size(); ++f)
1007 : : Assert( gid.size() == Lp[f].size(), "Size mismatch" );
1008 : :
1009 [ + + ]: 4236 : for (std::size_t i=0; i<gid.size(); ++i) {
1010 : : auto& nfu = m_uNodefieldsc[ gid[i] ];
1011 : 4164 : nfu.first.resize( Lu.size() );
1012 [ + + ]: 41640 : for (std::size_t f=0; f<Lu.size(); ++f) nfu.first[f] += Lu[f][i];
1013 : 4164 : nfu.second += nesup[i];
1014 : 4164 : auto& nfp = m_pNodefieldsc[ gid[i] ];
1015 : 4164 : nfp.first.resize( Lp.size() );
1016 [ + + ]: 24984 : for (std::size_t f=0; f<Lp.size(); ++f) nfp.first[f] += Lp[f][i];
1017 : 4164 : nfp.second += nesup[i];
1018 : : }
1019 : :
1020 : : // When we have heard from all chares we communicate with, this chare is done
1021 [ + + ]: 72 : if (++m_nnod == Disc()->NodeCommMap().size()) {
1022 : 24 : m_nnod = 0;
1023 : 24 : comnodeout_complete();
1024 : : }
1025 : 72 : }
1026 : :
1027 : : void
1028 : 1668 : FV::stage()
1029 : : // *****************************************************************************
1030 : : // Evaluate whether to continue with next time step stage
1031 : : // *****************************************************************************
1032 : : {
1033 : : // Increment Runge-Kutta stage counter
1034 : 1668 : ++m_stage;
1035 : :
1036 : : // if not all Runge-Kutta stages complete, continue to next time stage,
1037 : : // otherwise prepare for nodal field output
1038 [ + + ]: 1668 : if (m_stage < m_nrk)
1039 : 834 : next();
1040 : : else
1041 [ + - ][ + - ]: 2502 : startFieldOutput( CkCallback(CkIndex_FV::step(), thisProxy[thisIndex]) );
[ - + ][ - - ]
1042 : 1668 : }
1043 : :
1044 : : void
1045 : 536 : FV::evalLB( int nrestart )
1046 : : // *****************************************************************************
1047 : : // Evaluate whether to do load balancing
1048 : : //! \param[in] nrestart Number of times restarted
1049 : : // *****************************************************************************
1050 : : {
1051 : 536 : auto d = Disc();
1052 : :
1053 : : // Detect if just returned from a checkpoint and if so, zero timers and flag
1054 [ - + ]: 536 : if (d->restarted( nrestart )) m_finished = 0;
1055 : :
1056 : 536 : const auto lbfreq = g_inputdeck.get< tag::cmd, tag::lbfreq >();
1057 : 536 : const auto nonblocking = g_inputdeck.get< tag::cmd, tag::nonblocking >();
1058 : :
1059 : : // Load balancing if user frequency is reached or after the second time-step
1060 [ - + ][ - - ]: 536 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
1061 : :
1062 : 536 : AtSync();
1063 [ - + ]: 536 : if (nonblocking) next();
1064 : :
1065 : : } else {
1066 : :
1067 : 0 : next();
1068 : :
1069 : : }
1070 : 536 : }
1071 : :
1072 : : void
1073 : 536 : FV::evalRestart()
1074 : : // *****************************************************************************
1075 : : // Evaluate whether to save checkpoint/restart
1076 : : // *****************************************************************************
1077 : : {
1078 : 536 : auto d = Disc();
1079 : :
1080 : 536 : const auto rsfreq = g_inputdeck.get< tag::cmd, tag::rsfreq >();
1081 : 536 : const auto benchmark = g_inputdeck.get< tag::cmd, tag::benchmark >();
1082 : :
1083 [ + + ][ - + ]: 536 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
1084 : :
1085 : 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
1086 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
1087 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
[ - - ]
1088 : :
1089 : : } else {
1090 : :
1091 : 536 : evalLB( /* nrestart = */ -1 );
1092 : :
1093 : : }
1094 : 536 : }
1095 : :
1096 : : void
1097 : 834 : FV::step()
1098 : : // *****************************************************************************
1099 : : // Evaluate wether to continue with next time step
1100 : : // *****************************************************************************
1101 : : {
1102 : 834 : auto d = Disc();
1103 : :
1104 : : // Output time history
1105 [ + - ][ + - ]: 834 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ]
1106 : 0 : std::vector< std::vector< tk::real > > hist;
1107 [ - - ][ - - ]: 0 : auto h = g_fvpde[d->MeshId()].histOutput( d->Hist(), myGhosts()->m_inpoel,
1108 [ - - ][ - - ]: 0 : myGhosts()->m_coord, m_u, m_p );
1109 [ - - ]: 0 : hist.insert( end(hist), begin(h), end(h) );
1110 [ - - ]: 0 : d->history( std::move(hist) );
1111 : : }
1112 : :
1113 : : // Output one-liner status report to screen
1114 : 834 : d->status();
1115 : : // Reset Runge-Kutta stage counter
1116 : 834 : m_stage = 0;
1117 : :
1118 : : // If neither max iterations nor max time reached, continue, otherwise finish
1119 [ + + ]: 834 : if (not m_finished) {
1120 : :
1121 : 536 : evalRestart();
1122 : :
1123 : : } else {
1124 : :
1125 : 298 : auto meshid = d->MeshId();
1126 [ + - ]: 596 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1127 : 596 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
1128 : :
1129 : : }
1130 : 834 : }
1131 : :
1132 : : #include "NoWarning/fv.def.h"
|
