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