Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/DG.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief DG advances a system of PDEs with the discontinuous Galerkin scheme
9 : : \details DG advances a system of partial differential equations (PDEs) using
10 : : discontinuous Galerkin (DG) finite element (FE) spatial discretization (on
11 : : tetrahedron elements) combined with Runge-Kutta (RK) time stepping.
12 : : \see The documentation in DG.h.
13 : : */
14 : : // *****************************************************************************
15 : :
16 : : #include <algorithm>
17 : : #include <numeric>
18 : : #include <sstream>
19 : :
20 : : #include "DG.hpp"
21 : : #include "Discretization.hpp"
22 : : #include "DGPDE.hpp"
23 : : #include "DiagReducer.hpp"
24 : : #include "DerivedData.hpp"
25 : : #include "ElemDiagnostics.hpp"
26 : : #include "Inciter/InputDeck/InputDeck.hpp"
27 : : #include "Refiner.hpp"
28 : : #include "Limiter.hpp"
29 : : #include "Reorder.hpp"
30 : : #include "Vector.hpp"
31 : : #include "Around.hpp"
32 : : #include "Integrate/Basis.hpp"
33 : : #include "FieldOutput.hpp"
34 : : #include "ChareStateCollector.hpp"
35 : :
36 : : namespace inciter {
37 : :
38 : : extern ctr::InputDeck g_inputdeck;
39 : : extern ctr::InputDeck g_inputdeck_defaults;
40 : : extern std::vector< DGPDE > g_dgpde;
41 : :
42 : : //! Runge-Kutta coefficients
43 : : static const std::array< std::array< tk::real, 3 >, 2 >
44 : : rkcoef{{ {{ 0.0, 3.0/4.0, 1.0/3.0 }}, {{ 1.0, 1.0/4.0, 2.0/3.0 }} }};
45 : :
46 : : } // inciter::
47 : :
48 : : extern tk::CProxy_ChareStateCollector stateProxy;
49 : :
50 : : using inciter::DG;
51 : :
52 : 918 : DG::DG( const CProxy_Discretization& disc,
53 : : const std::map< int, std::vector< std::size_t > >& bface,
54 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
55 : 918 : const std::vector< std::size_t >& triinpoel ) :
56 : : m_disc( disc ),
57 : : m_ndof_NodalExtrm( 3 ), // for the first order derivatives in 3 directions
58 : : m_ncomfac( 0 ),
59 : : m_nadj( 0 ),
60 : : m_ncomEsup( 0 ),
61 : : m_nsol( 0 ),
62 : : m_ninitsol( 0 ),
63 : : m_nlim( 0 ),
64 : : m_nnod( 0 ),
65 : : m_nreco( 0 ),
66 : : m_nnodalExtrema( 0 ),
67 : : m_inpoel( Disc()->Inpoel() ),
68 : 918 : m_coord( Disc()->Coord() ),
69 [ + - ]: 1836 : m_fd( m_inpoel, bface, tk::remap(triinpoel,Disc()->Lid()) ),
70 : 918 : m_u( m_inpoel.size()/4,
71 : 1836 : g_inputdeck.get< tag::discr, tag::rdof >()*
72 : 918 : g_inputdeck.get< tag::component >().nprop() ),
73 : : m_un( m_u.nunk(), m_u.nprop() ),
74 : : m_p( m_u.nunk(),
75 : 1836 : g_inputdeck.get< tag::discr, tag::rdof >()*
76 [ + - ]: 918 : std::accumulate( begin(g_dgpde), end(g_dgpde), 0u,
77 : 918 : [](std::size_t s, const DGPDE& eq){ return s + eq.nprim(); } ) ),
78 : 918 : m_geoFace( tk::genGeoFaceTri( m_fd.Nipfac(), m_fd.Inpofa(), m_coord) ),
79 : 918 : m_geoElem( tk::genGeoElemTet( m_inpoel, m_coord ) ),
80 : : m_lhs( m_u.nunk(),
81 : 1836 : g_inputdeck.get< tag::discr, tag::ndof >()*
82 : 918 : g_inputdeck.get< tag::component >().nprop() ),
83 : : m_rhs( m_u.nunk(), m_lhs.nprop() ),
84 : : m_uNodalExtrm(),
85 : : m_pNodalExtrm(),
86 : : m_uNodalExtrmc(),
87 : : m_pNodalExtrmc(),
88 : 918 : m_nfac( m_fd.Inpofa().size()/3 ),
89 : 918 : m_nunk( m_u.nunk() ),
90 : 918 : m_npoin( m_coord[0].size() ),
91 : : m_ipface(),
92 : : m_bndFace(),
93 : : m_ghostData(),
94 : : m_sendGhost(),
95 : : m_ghostReq( 0 ),
96 : : m_ghost(),
97 : : m_exptGhost(),
98 : : m_recvGhost(),
99 : : m_diag(),
100 : : m_stage( 0 ),
101 : : m_ndof(),
102 : : m_numEqDof(),
103 : : m_bid(),
104 : : m_uc(),
105 : : m_pc(),
106 : : m_ndofc(),
107 : : m_initial( 1 ),
108 : : m_expChBndFace(),
109 : : m_infaces(),
110 : : m_esup(),
111 : : m_esupc(),
112 : : m_elemfields(),
113 : : m_nodefields(),
114 : : m_nodefieldsc(),
115 : : m_outmesh(),
116 : : m_boxelems(),
117 [ + - ][ + - ]: 9180 : m_shockmarker(m_u.nunk())
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
118 : : // *****************************************************************************
119 : : // Constructor
120 : : //! \param[in] disc Discretization proxy
121 : : //! \param[in] bface Boundary-faces mapped to side set ids
122 : : //! \param[in] triinpoel Boundary-face connectivity
123 : : // *****************************************************************************
124 : : {
125 [ + + ][ + + ]: 1798 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
126 [ + + ]: 880 : g_inputdeck.get< tag::cmd, tag::quiescence >())
127 [ + - ][ + - ]: 541 : stateProxy.ckLocalBranch()->insert( "DG", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
[ + - ]
128 : : "DG" );
129 : :
130 : : // assign number of dofs for each equation in all pde systems
131 [ + + ]: 1836 : for (const auto& eq : g_dgpde) {
132 [ + - ]: 918 : eq.numEquationDofs(m_numEqDof);
133 : : }
134 : :
135 : : // Allocate storage for the vector of nodal extrema
136 [ + - ][ + - ]: 1836 : m_uNodalExtrm.resize( Disc()->Bid().size(), std::vector<tk::real>( 2*
137 [ + - ]: 918 : m_ndof_NodalExtrm*g_inputdeck.get< tag::component >().nprop() ) );
138 [ + - ][ + - ]: 1836 : m_pNodalExtrm.resize( Disc()->Bid().size(), std::vector<tk::real>( 2*
139 [ + - ]: 918 : m_ndof_NodalExtrm*m_p.nprop()/g_inputdeck.get< tag::discr, tag::rdof >()));
140 : :
141 : : // Initialization for the buffer vector of nodal extrema
142 [ + - ]: 918 : resizeNodalExtremac();
143 : :
144 : 918 : usesAtSync = true; // enable migration at AtSync
145 : :
146 : : // Enable SDAG wait for setting up chare boundary faces
147 [ + - ][ + - ]: 918 : thisProxy[ thisIndex ].wait4fac();
148 : :
149 : : // Ensure that mesh partition is not leaky
150 [ + - ][ - + ]: 918 : Assert( !tk::leakyPartition(m_fd.Esuel(), m_inpoel, m_coord),
[ - - ][ - - ]
[ - - ]
151 : : "Input mesh to DG leaky" );
152 : :
153 : : // Ensure mesh physical boundary for the entire problem not leaky,
154 : : // effectively checking if the user has specified boundary conditions on all
155 : : // physical boundary faces
156 [ + - ]: 918 : bndIntegral();
157 : 918 : }
158 : :
159 : : void
160 : 918 : DG::bndIntegral()
161 : : // *****************************************************************************
162 : : // Compute partial boundary surface integral and sum across all chares
163 : : //! \details This function computes a partial surface integral over the boundary
164 : : //! of the faces of this mesh partition then sends its contribution to perform
165 : : //! the integral acorss the total problem boundary. After the global sum a
166 : : //! non-zero vector result indicates a leak, e.g., a hole in the boundary
167 : : //! which indicates an error in the boundary face data structures used to
168 : : //! compute the partial surface integrals.
169 : : // *****************************************************************************
170 : : {
171 : : // Storage for surface integral over our mesh chunk physical boundary
172 [ + - ]: 918 : std::vector< tk::real > s{{ 0.0, 0.0, 0.0 }};
173 : :
174 : : // Integrate over all physical boundary faces
175 [ + - ][ + + ]: 232884 : for (std::size_t f=0; f<m_fd.Nbfac(); ++f) {
176 [ + - ][ + - ]: 231966 : s[0] += m_geoFace(f,0,0) * m_geoFace(f,1,0);
177 [ + - ][ + - ]: 231966 : s[1] += m_geoFace(f,0,0) * m_geoFace(f,2,0);
178 [ + - ][ + - ]: 231966 : s[2] += m_geoFace(f,0,0) * m_geoFace(f,3,0);
179 : : }
180 : :
181 [ + - ]: 918 : s.push_back( 1.0 ); // positive: call-back to resizeComm() after reduction
182 [ + - ][ + - ]: 918 : s.push_back( static_cast< tk::real >( Disc()->MeshId() ) );
183 : :
184 : : // Send contribution to host summing partial surface integrals
185 [ + - ]: 918 : contribute( s, CkReduction::sum_double,
186 [ + - ][ + - ]: 1836 : CkCallback(CkReductionTarget(Transporter,bndint), Disc()->Tr()) );
[ + - ]
187 : 918 : }
188 : :
189 : : void
190 : 1026 : DG::resizeComm()
191 : : // *****************************************************************************
192 : : // Start sizing communication buffers and setting up ghost data
193 : : // *****************************************************************************
194 : : {
195 [ + - ]: 1026 : auto d = Disc();
196 : :
197 : 1026 : const auto& gid = d->Gid();
198 : 1026 : const auto& inpofa = m_fd.Inpofa();
199 : 1026 : const auto& esuel = m_fd.Esuel();
200 : :
201 : : // Perform leak test on mesh partition
202 [ + - ][ - + ]: 1026 : Assert( !tk::leakyPartition( esuel, m_inpoel, m_coord ),
[ - - ][ - - ]
[ - - ]
203 : : "Mesh partition leaky" );
204 : :
205 : : // Activate SDAG waits for face adjacency map (ghost data) calculation
206 [ + - ][ + - ]: 1026 : thisProxy[ thisIndex ].wait4ghost();
207 [ + - ][ + - ]: 1026 : thisProxy[ thisIndex ].wait4esup();
208 : :
209 : : // Enable SDAG wait for initially building the solution vector and limiting
210 [ + + ]: 1026 : if (m_initial) {
211 [ + - ][ + - ]: 918 : thisProxy[ thisIndex ].wait4sol();
212 [ + - ][ + - ]: 918 : thisProxy[ thisIndex ].wait4reco();
213 [ + - ][ + - ]: 918 : thisProxy[ thisIndex ].wait4nodalExtrema();
214 [ + - ][ + - ]: 918 : thisProxy[ thisIndex ].wait4lim();
215 [ + - ][ + - ]: 918 : thisProxy[ thisIndex ].wait4nod();
216 : : }
217 : :
218 : : // Invert inpofa to enable searching for faces based on (global) node triplets
219 [ - + ][ - - ]: 1026 : Assert( inpofa.size() % 3 == 0, "Inpofa must contain triplets" );
[ - - ][ - - ]
220 [ + + ]: 1837956 : for (std::size_t f=0; f<inpofa.size()/3; ++f)
221 [ + - ]: 3673860 : m_ipface.insert( {{{ gid[ inpofa[f*3+0] ],
222 : 1836930 : gid[ inpofa[f*3+1] ],
223 : 1836930 : gid[ inpofa[f*3+2] ] }}} );
224 : :
225 : : // At this point ipface has node-id-triplets (faces) on the internal
226 : : // chare-domain and on the physical boundary but not on chare boundaries,
227 : : // hence the name internal + physical boundary faces.
228 : :
229 : : // Build a set of faces (each face given by 3 global node IDs) associated to
230 : : // chares we potentially share boundary faces with.
231 : 2052 : tk::UnsMesh::FaceSet potbndface;
232 [ + + ]: 892757 : for (std::size_t e=0; e<esuel.size()/4; ++e) { // for all our tets
233 : 891731 : auto mark = e*4;
234 [ + + ]: 4458655 : for (std::size_t f=0; f<4; ++f) // for all tet faces
235 [ + + ]: 3566924 : if (esuel[mark+f] == -1) { // if face has no outside-neighbor tet
236 : : // if does not exist among the internal and physical boundary faces,
237 : : // store as a potential chare-boundary face
238 : 379396 : tk::UnsMesh::Face t{{ gid[ m_inpoel[ mark + tk::lpofa[f][0] ] ],
239 : 379396 : gid[ m_inpoel[ mark + tk::lpofa[f][1] ] ],
240 : 758792 : gid[ m_inpoel[ mark + tk::lpofa[f][2] ] ] }};
241 [ + - ][ + + ]: 379396 : if (m_ipface.find(t) == end(m_ipface)) {
242 [ + - ][ - + ]: 136230 : Assert( m_expChBndFace.insert(t).second,
[ - - ][ - - ]
[ - - ]
243 : : "Store expected chare-boundary face" );
244 [ + - ]: 136230 : potbndface.insert( t );
245 : : }
246 : : }
247 : : }
248 : :
249 [ - + ][ - - ]: 1026 : if ( g_inputdeck.get< tag::cmd, tag::feedback >() ) d->Tr().chbndface();
250 : :
251 : : // In the following we assume that the size of the (potential) boundary-face
252 : : // adjacency map above does not necessarily equal to that of the node
253 : : // adjacency map. This is because while a node can be shared at a single
254 : : // corner or along an edge, that does not necessarily share a face as well
255 : : // (in other words, shared nodes or edges can exist that are not part of a
256 : : // shared face). So the chares we communicate with across faces are not
257 : : // necessarily the same as the chares we would communicate nodes with.
258 : : //
259 : : // Since the sizes of the node and face adjacency maps are not the same, while
260 : : // sending the faces on chare boundaries would be okay, however, the receiver
261 : : // would not necessarily know how many chares it must receive from. To solve
262 : : // this problem we send to chares which we share at least a single node with,
263 : : // i.e., rely on the node-adjacency map. Note that to all chares we share at
264 : : // least a single node with we send all our potential chare-boundary faces.
265 : : // This is the same list of faces to all chares we send.
266 : : //
267 : : // Another underlying assumption here is, of course, that the size of the face
268 : : // adjacency map is always smaller than or equal to that of the node adjacency
269 : : // map, which is always true. Since the receive side already knows how many
270 : : // fellow chares it must receive shared node ids from, we use that to detect
271 : : // completion of the number of receives in comfac(). This simplifies the
272 : : // communication pattern and code.
273 : :
274 : : // Send sets of faces adjacent to chare boundaries to fellow workers (if any)
275 [ + + ]: 1026 : if (d->NodeCommMap().empty()) // in serial, skip setting up ghosts altogether
276 [ + - ]: 43 : faceAdj();
277 : : else
278 : : // for all chares we share nodes with
279 [ + + ]: 9211 : for (const auto& c : d->NodeCommMap()) {
280 [ + - ][ + - ]: 8228 : thisProxy[ c.first ].comfac( thisIndex, potbndface );
281 : : }
282 : :
283 [ + - ]: 1026 : ownfac_complete();
284 : 1026 : }
285 : :
286 : : bool
287 : 1026 : DG::leakyAdjacency()
288 : : // *****************************************************************************
289 : : // Perform leak-test on chare boundary faces
290 : : //! \details This function computes a surface integral over the boundary of the
291 : : //! faces after the face adjacency communication map is completed. A non-zero
292 : : //! vector result indicates a leak, e.g., a hole in the partition (covered by
293 : : //! the faces of the face adjacency communication map), which indicates an
294 : : //! error upstream in the code that sets up the face communication data
295 : : //! structures.
296 : : //! \note Compared to tk::leakyPartition() this function performs the leak-test
297 : : //! on the face geometry data structure enlarged by ghost faces on this
298 : : //! partition by computing a discrete surface integral considering the
299 : : //! physical and chare boundary faces, which should be equal to zero for a
300 : : //! closed domain.
301 : : //! \return True if our chare face adjacency leaks.
302 : : // *****************************************************************************
303 : : {
304 : : // Storage for surface integral over our chunk of the adjacency
305 : 1026 : std::array< tk::real, 3 > s{{ 0.0, 0.0, 0.0 }};
306 : :
307 : : // physical boundary faces
308 [ + - ][ + + ]: 244192 : for (std::size_t f=0; f<m_fd.Nbfac(); ++f) {
309 [ + - ][ + - ]: 243166 : s[0] += m_geoFace(f,0,0) * m_geoFace(f,1,0);
310 [ + - ][ + - ]: 243166 : s[1] += m_geoFace(f,0,0) * m_geoFace(f,2,0);
311 [ + - ][ + - ]: 243166 : s[2] += m_geoFace(f,0,0) * m_geoFace(f,3,0);
312 : : }
313 : :
314 : : // chare-boundary faces
315 [ + + ]: 137256 : for (std::size_t f=m_fd.Nipfac(); f<m_fd.Esuf().size()/2; ++f) {
316 [ + - ][ + - ]: 136230 : s[0] += m_geoFace(f,0,0) * m_geoFace(f,1,0);
317 [ + - ][ + - ]: 136230 : s[1] += m_geoFace(f,0,0) * m_geoFace(f,2,0);
318 [ + - ][ + - ]: 136230 : s[2] += m_geoFace(f,0,0) * m_geoFace(f,3,0);
319 : : }
320 : :
321 : 1026 : auto eps = std::numeric_limits< tk::real >::epsilon() * 100;
322 [ + - ][ + - ]: 1026 : return std::abs(s[0]) > eps || std::abs(s[1]) > eps || std::abs(s[2]) > eps;
[ - + ]
323 : : }
324 : :
325 : : bool
326 : 1026 : DG::faceMatch()
327 : : // *****************************************************************************
328 : : // Check if esuf of chare-boundary faces matches
329 : : //! \details This function checks each chare-boundary esuf entry for the left
330 : : //! and right elements. Then, it tries to match all vertices of these
331 : : //! elements. Exactly three of these vertices must match if the esuf entry
332 : : //! has been updated correctly at chare-boundaries.
333 : : //! \return True if chare-boundary faces match.
334 : : // *****************************************************************************
335 : : {
336 : 1026 : const auto& esuf = m_fd.Esuf();
337 : 1026 : bool match(true);
338 : :
339 : 1026 : auto eps = std::numeric_limits< tk::real >::epsilon() * 100;
340 : :
341 [ + + ]: 137256 : for (auto f=m_fd.Nipfac(); f<esuf.size()/2; ++f)
342 : : {
343 : 136230 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
344 : 136230 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
345 : :
346 : 136230 : std::size_t count = 0;
347 : :
348 [ + + ]: 681150 : for (std::size_t i=0; i<4; ++i)
349 : : {
350 : 544920 : auto ip = m_inpoel[4*el+i];
351 [ + + ]: 2724600 : for (std::size_t j=0; j<4; ++j)
352 : : {
353 : 2179680 : auto jp = m_inpoel[4*er+j];
354 : 2179680 : auto xdiff = std::abs( m_coord[0][ip] - m_coord[0][jp] );
355 : 2179680 : auto ydiff = std::abs( m_coord[1][ip] - m_coord[1][jp] );
356 : 2179680 : auto zdiff = std::abs( m_coord[2][ip] - m_coord[2][jp] );
357 : :
358 [ + + ][ + + ]: 2179680 : if ( xdiff<=eps && ydiff<=eps && zdiff<=eps ) ++count;
[ + + ]
359 : : }
360 : : }
361 : :
362 [ + - ][ + - ]: 136230 : match = (match && count == 3);
363 : : }
364 : :
365 : 1026 : return match;
366 : : }
367 : :
368 : : void
369 : 8228 : DG::comfac( int fromch, const tk::UnsMesh::FaceSet& infaces )
370 : : // *****************************************************************************
371 : : // Receive unique set of faces we potentially share with/from another chare
372 : : //! \param[in] fromch Sender chare id
373 : : //! \param[in] infaces Unique set of faces we potentially share with fromch
374 : : // *****************************************************************************
375 : : {
376 [ + + ][ + + ]: 15990 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
377 [ + + ]: 7762 : g_inputdeck.get< tag::cmd, tag::quiescence >())
378 [ + - ][ + - ]: 4840 : stateProxy.ckLocalBranch()->insert( "DG", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
379 : : "comfac" );
380 : :
381 : : // Buffer up incoming data
382 : 8228 : m_infaces[ fromch ] = infaces;
383 : :
384 : : // if we have heard from all fellow chares that we share at least a single
385 : : // node, edge, or face with
386 [ + + ]: 8228 : if (++m_ncomfac == Disc()->NodeCommMap().size()) {
387 : 983 : m_ncomfac = 0;
388 : 983 : comfac_complete();
389 : : }
390 : 8228 : }
391 : :
392 : : void
393 : 983 : DG::bndFaces()
394 : : // *****************************************************************************
395 : : // Compute chare-boundary faces
396 : : //! \details This is called when both send and receives are completed on a
397 : : //! chare and thus we are ready to compute chare-boundary faces and ghost data.
398 : : // *****************************************************************************
399 : : {
400 : 983 : auto d = Disc();
401 [ - + ]: 983 : if ( g_inputdeck.get< tag::cmd, tag::feedback >() ) d->Tr().chcomfac();
402 : 983 : const auto& esuel = m_fd.Esuel();
403 : 983 : const auto& gid = d->Gid();
404 : :
405 [ + + ]: 9211 : for (const auto& in : m_infaces) {
406 : : // Find sender chare among chares we potentially share faces with. Note that
407 : : // it is feasible that a sender chare called us but we do not have a set of
408 : : // faces associated to that chare. This can happen if we only share a single
409 : : // node or an edge but not a face with that chare.
410 [ + - ]: 8228 : auto& bndface = m_bndFace[ in.first ]; // will associate to sender chare
411 : : // Try to find incoming faces on our chare boundary with other chares. If
412 : : // found, generate and assign new local face ID, associated to sender chare.
413 [ + + ]: 3886374 : for (std::size_t e=0; e<esuel.size()/4; ++e) { // for all our tets
414 : 3878146 : auto mark = e*4;
415 [ + + ]: 19390730 : for (std::size_t f=0; f<4; ++f) { // for all cell faces
416 [ + + ]: 15512584 : if (esuel[mark+f] == -1) { // if face has no outside-neighbor tet
417 : 2014270 : tk::UnsMesh::Face t{{ gid[ m_inpoel[ mark + tk::lpofa[f][0] ] ],
418 : 2014270 : gid[ m_inpoel[ mark + tk::lpofa[f][1] ] ],
419 : 4028540 : gid[ m_inpoel[ mark + tk::lpofa[f][2] ] ] }};
420 : : // if found among the incoming faces and if not one of our internal
421 : : // nor physical boundary faces
422 [ + - ][ + + ]: 2150500 : if ( in.second.find(t) != end(in.second) &&
[ + - ]
423 [ + - ][ + + ]: 2150500 : m_ipface.find(t) == end(m_ipface) ) {
424 [ + - ]: 136230 : bndface[t][0] = m_nfac++; // assign new local face ID
425 : : }
426 : : }
427 : : }
428 : : }
429 : : // If at this point if we have not found any face among our faces we
430 : : // potentially share with fromch, there is no need to keep an empty set of
431 : : // faces associated to fromch as we only share nodes or edges with it, but
432 : : // not faces.
433 [ + + ][ + - ]: 8228 : if (bndface.empty()) m_bndFace.erase( in.first );
434 : : }
435 : :
436 : 983 : tk::destroy(m_ipface);
437 : 983 : tk::destroy(m_infaces);
438 : :
439 : : // Ensure all expected faces have been received
440 [ - + ][ - - ]: 983 : Assert( receivedChBndFaces(),
[ - - ][ - - ]
441 : : "Expected and received chare boundary faces mismatch" );
442 : :
443 : : // Basic error checking on chare-boundary-face map
444 [ + - ][ - + ]: 983 : Assert( m_bndFace.find( thisIndex ) == m_bndFace.cend(),
[ - - ][ - - ]
[ - - ]
445 : : "Face-communication map should not contain data for own chare ID" );
446 : :
447 : : // Store (local) tet ID adjacent to our chare boundary from the inside
448 [ + + ]: 683566 : for (std::size_t e=0; e<esuel.size()/4; ++e) { // for all our tets
449 : 682583 : auto mark = e*4;
450 [ + + ]: 3412915 : for (std::size_t f=0; f<4; ++f) { // for all cell faces
451 [ + + ]: 2730332 : if (esuel[mark+f] == -1) { // if face has no outside-neighbor tet
452 : 309050 : tk::UnsMesh::Face t{{ gid[ m_inpoel[ mark + tk::lpofa[f][0] ] ],
453 : 309050 : gid[ m_inpoel[ mark + tk::lpofa[f][1] ] ],
454 : 618100 : gid[ m_inpoel[ mark + tk::lpofa[f][2] ] ] }};
455 [ + - ]: 309050 : auto c = findchare( t );
456 [ + + ]: 309050 : if (c > -1) {
457 [ + - ]: 136230 : auto& lbndface = tk::ref_find( m_bndFace, c );
458 [ + - ]: 136230 : auto& face = tk::ref_find( lbndface, t );
459 : 136230 : face[1] = e; // store (local) inner tet ID adjacent to face
460 : : }
461 : : }
462 : : }
463 : : }
464 : :
465 : : // At this point m_bndFace is complete on this PE. This means that starting
466 : : // from the sets of faces we potentially share with fellow chares we now
467 : : // only have those faces we actually share faces with (through which we need
468 : : // to communicate later). Also, m_bndFace not only has the unique faces
469 : : // associated to fellow chares, but also a newly assigned local face ID as
470 : : // well as the local id of the inner tet adjacent to the face. Continue by
471 : : // starting setting up ghost data
472 : 983 : setupGhost();
473 : : // Besides setting up our own ghost data, we also issue requests (for ghost
474 : : // data) to those chares which we share faces with. Note that similar to
475 : : // comfac() we are calling reqGhost() by going through the node communication
476 : : // map instead, which may send requests to those chare we do not share faces
477 : : // with. This is so that we can test for completing by querying the size of
478 : : // the already complete node commincation map in reqGhost. Requests in
479 : : // sendGhost will only be fullfilled based on m_ghostData.
480 [ + + ]: 9211 : for (const auto& c : d->NodeCommMap()) // for all chares we share nodes with
481 [ + - ][ + - ]: 8228 : thisProxy[ c.first ].reqGhost();
482 : 983 : }
483 : :
484 : : bool
485 : 983 : DG::receivedChBndFaces()
486 : : // *****************************************************************************
487 : : // Verify that all chare-boundary faces have been received
488 : : //! \return True if all chare-boundary faces have been received
489 : : // *****************************************************************************
490 : : {
491 [ + - ]: 983 : auto d = Disc();
492 : 1966 : tk::UnsMesh::FaceSet recvBndFace;
493 : :
494 : : // Collect chare-boundary faces that have been received and expected
495 [ + + ]: 6165 : for (const auto& c : m_bndFace)
496 [ + + ]: 141412 : for (const auto& f : c.second)
497 [ + - ][ + - ]: 136230 : if (m_expChBndFace.find(f.first) != end(m_expChBndFace))
498 [ + - ]: 136230 : recvBndFace.insert(f.first);
499 : :
500 : : // Collect info on expected but not received faces
501 [ + - ]: 1966 : std::stringstream msg;
502 [ + + ]: 137213 : for (const auto& f : m_expChBndFace)
503 [ + - ][ - + ]: 136230 : if (recvBndFace.find(f) == end(recvBndFace)) {
504 : 0 : const auto& x = m_coord[0];
505 : 0 : const auto& y = m_coord[1];
506 : 0 : const auto& z = m_coord[2];
507 [ - - ]: 0 : auto A = tk::cref_find( d->Lid(), f[0] );
508 [ - - ]: 0 : auto B = tk::cref_find( d->Lid(), f[1] );
509 [ - - ]: 0 : auto C = tk::cref_find( d->Lid(), f[2] );
510 [ - - ][ - - ]: 0 : msg << '{' << A << ',' << B << ',' << C << "}:("
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
511 [ - - ][ - - ]: 0 : << x[A] << ',' << y[A] << ',' << z[A] << ' '
[ - - ][ - - ]
[ - - ][ - - ]
512 [ - - ][ - - ]: 0 : << x[B] << ',' << y[B] << ',' << z[B] << ' '
[ - - ][ - - ]
[ - - ][ - - ]
513 [ - - ][ - - ]: 0 : << x[C] << ',' << y[C] << ',' << z[C] << ") ";
[ - - ][ - - ]
[ - - ][ - - ]
514 : : }
515 : :
516 : 983 : tk::destroy( m_expChBndFace );
517 : :
518 : : // Error out with info on missing faces
519 [ + - ]: 983 : auto s = msg.str();
520 [ - + ]: 983 : if (!s.empty()) {
521 [ - - ][ - - ]: 0 : Throw( "DG chare " + std::to_string(thisIndex) +
[ - - ][ - - ]
[ - - ][ - - ]
522 : : " missing face(s) {local node ids} (node coords): " + s );
523 : : } else {
524 : 1966 : return true;
525 : : }
526 : : }
527 : :
528 : : void
529 : 983 : DG::setupGhost()
530 : : // *****************************************************************************
531 : : // Setup own ghost data on this chare
532 : : // *****************************************************************************
533 : : {
534 : 983 : auto d = Disc();
535 : 983 : const auto& gid = d->Gid();
536 : :
537 : : // Enlarge elements surrounding faces data structure for ghosts
538 [ + - ]: 983 : m_fd.Esuf().resize( 2*m_nfac, -2 );
539 [ + - ]: 983 : m_fd.Inpofa().resize( 3*m_nfac, 0 );
540 : : // Enlarge face geometry data structure for ghosts
541 : 983 : m_geoFace.resize( m_nfac, 0.0 );
542 : :
543 : 983 : const auto& esuel = m_fd.Esuel();
544 : :
545 : : // Collect tet ids, their face connectivity (given by 3 global node IDs, each
546 : : // triplet for potentially multiple faces on the chare boundary), and their
547 : : // elem geometry data (see GhostData) associated to fellow chares adjacent to
548 : : // chare boundaries. Once received by fellow chares, these tets will become
549 : : // known as ghost elements and their data as ghost data.
550 [ + + ]: 683566 : for (std::size_t e=0; e<esuel.size()/4; ++e) { // for all our tets
551 : 682583 : auto mark = e*4;
552 [ + + ]: 3412915 : for (std::size_t f=0; f<4; ++f) { // for all cell faces
553 [ + + ]: 2730332 : if (esuel[mark+f] == -1) { // if face has no outside-neighbor tet
554 : 309050 : tk::UnsMesh::Face t{{ gid[ m_inpoel[ mark + tk::lpofa[f][0] ] ],
555 : 309050 : gid[ m_inpoel[ mark + tk::lpofa[f][1] ] ],
556 : 618100 : gid[ m_inpoel[ mark + tk::lpofa[f][2] ] ] }};
557 [ + - ]: 309050 : auto c = findchare( t );
558 : : // It is possible that we do not find the chare for this face. We are
559 : : // looping through all of our tets and interrogating all faces that do
560 : : // not have neighboring tets but we only care about chare-boundary faces
561 : : // here as only those need ghost data. (esuel may also contain
562 : : // physical boundary faces)
563 [ + + ]: 309050 : if (c > -1) {
564 : : // Will store ghost data associated to neighbor chare
565 [ + - ]: 136230 : auto& ghost = m_ghostData[ c ];
566 : : // Store tet id adjacent to chare boundary as key for ghost data
567 [ + - ]: 136230 : auto& tuple = ghost[ e ];
568 : : // If tetid e has not yet been encountered, store geometry (only once)
569 : 136230 : auto& nodes = std::get< 0 >( tuple );
570 [ + + ]: 136230 : if (nodes.empty()) {
571 [ + - ]: 106332 : std::get< 1 >( tuple ) = m_geoElem[ e ];
572 : :
573 : 106332 : auto& ncoord = std::get< 2 >( tuple );
574 : 106332 : ncoord[0] = m_coord[0][ m_inpoel[ mark+f ] ];
575 : 106332 : ncoord[1] = m_coord[1][ m_inpoel[ mark+f ] ];
576 : 106332 : ncoord[2] = m_coord[2][ m_inpoel[ mark+f ] ];
577 : :
578 : 106332 : std::get< 3 >( tuple ) = f;
579 : :
580 : 212664 : std::get< 4 >( tuple ) = {{ gid[ m_inpoel[ mark ] ],
581 : 106332 : gid[ m_inpoel[ mark+1 ] ],
582 : 106332 : gid[ m_inpoel[ mark+2 ] ],
583 : 106332 : gid[ m_inpoel[ mark+3 ] ] }};
584 : : }
585 : : // (Always) store face node IDs on chare boundary, even if tetid e has
586 : : // already been stored. Thus we store potentially multiple faces along
587 : : // the same chare-boundary. This happens, e.g., when the boundary
588 : : // between chares is zig-zaggy enough to have 2 or even 3 faces of the
589 : : // same tet.
590 [ + - ]: 136230 : nodes.push_back( t[0] );
591 [ + - ]: 136230 : nodes.push_back( t[1] );
592 [ + - ]: 136230 : nodes.push_back( t[2] );
593 [ - + ][ - - ]: 136230 : Assert( nodes.size() <= 4*3, "Overflow of faces/tet to send" );
[ - - ][ - - ]
594 : : }
595 : : }
596 : : }
597 : : }
598 : :
599 : : // Basic error checking on local ghost data
600 [ + - ][ - + ]: 983 : Assert( m_ghostData.find( thisIndex ) == m_ghostData.cend(),
[ - - ][ - - ]
[ - - ]
601 : : "Chare-node adjacency map should not contain data for own chare ID" );
602 : :
603 : : // More in-depth error checking on local ghost data
604 [ + + ]: 6165 : for (const auto& c : m_ghostData)
605 [ + + ]: 111514 : for ([[maybe_unused]] const auto& t : c.second) {
606 [ - + ][ - - ]: 106332 : Assert( !std::get< 0 >( t.second ).empty(),
[ - - ][ - - ]
607 : : "Emtpy face vector in ghost data" );
608 [ - + ][ - - ]: 106332 : Assert( std::get< 0 >( t.second ).size() % 3 == 0,
[ - - ][ - - ]
609 : : "Face node IDs must be triplets" );
610 [ - + ][ - - ]: 106332 : Assert( std::get< 0 >( t.second ).size() <= 4*3, // <= 4*3 (4*numfaces)
[ - - ][ - - ]
611 : : "Max number of faces for a single ghost tet is 4" );
612 [ - + ][ - - ]: 106332 : Assert( !std::get< 1 >( t.second ).empty(),
[ - - ][ - - ]
613 : : "No elem geometry data for ghost" );
614 [ - + ][ - - ]: 106332 : Assert( std::get< 1 >( t.second ).size() == m_geoElem.nprop(),
[ - - ][ - - ]
615 : : "Elem geometry data for ghost must be for single tet" );
616 [ - + ][ - - ]: 106332 : Assert( !std::get< 2 >( t.second ).empty(),
[ - - ][ - - ]
617 : : "No nodal coordinate data for ghost" );
618 : : }
619 : :
620 : 983 : ownghost_complete();
621 : 983 : }
622 : :
623 : : void
624 : 8228 : DG::reqGhost()
625 : : // *****************************************************************************
626 : : // Receive requests for ghost data
627 : : // *****************************************************************************
628 : : {
629 [ + + ][ + + ]: 15990 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
630 [ + + ]: 7762 : g_inputdeck.get< tag::cmd, tag::quiescence >())
631 [ + - ][ + - ]: 4840 : stateProxy.ckLocalBranch()->insert( "DG", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
632 : : "reqGhost" );
633 : :
634 : : // If every chare we communicate with has requested ghost data from us, we may
635 : : // fulfill the requests, but only if we have already setup our ghost data.
636 [ + + ]: 8228 : if (++m_ghostReq == Disc()->NodeCommMap().size()) {
637 : 983 : m_ghostReq = 0;
638 : 983 : reqghost_complete();
639 : : }
640 : 8228 : }
641 : :
642 : : void
643 : 983 : DG::sendGhost()
644 : : // *****************************************************************************
645 : : // Send all of our ghost data to fellow chares
646 : : // *****************************************************************************
647 : : {
648 [ + + ][ + + ]: 1928 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
649 [ + + ]: 945 : g_inputdeck.get< tag::cmd, tag::quiescence >())
650 [ + - ][ + - ]: 555 : stateProxy.ckLocalBranch()->insert( "DG", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
651 : : "sendGhost" );
652 : :
653 [ + + ]: 6165 : for (const auto& c : m_ghostData)
654 [ + - ][ + - ]: 5182 : thisProxy[ c.first ].comGhost( thisIndex, c.second );
655 : :
656 [ - + ]: 983 : if ( g_inputdeck.get< tag::cmd, tag::feedback >() ) Disc()->Tr().chghost();
657 : 983 : }
658 : :
659 : : int
660 : 618100 : DG::findchare( const tk::UnsMesh::Face& t )
661 : : // *****************************************************************************
662 : : // Find any chare for face (given by 3 global node IDs)
663 : : //! \param[in] t Face given by three global node IDs
664 : : //! \return Chare ID if found among any of the chares we communicate along
665 : : //! faces with, -1 if the face cannot be found.
666 : : // *****************************************************************************
667 : : {
668 [ + + ]: 2714492 : for (const auto& cf : m_bndFace)
669 : : // cppcheck-suppress useStlAlgorithm
670 [ + - ][ + + ]: 2368852 : if (cf.second.find(t) != end(cf.second))
671 : 272460 : return cf.first;
672 : 345640 : return -1;
673 : : }
674 : :
675 : : void
676 : 5182 : DG::comGhost( int fromch, const GhostData& ghost )
677 : : // *****************************************************************************
678 : : // Receive ghost data on chare boundaries from fellow chare
679 : : //! \param[in] fromch Caller chare ID
680 : : //! \param[in] ghost Ghost data, see Inciter/FaceData.h for the type
681 : : // *****************************************************************************
682 : : {
683 [ + + ][ + + ]: 10124 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
684 [ + + ]: 4942 : g_inputdeck.get< tag::cmd, tag::quiescence >())
685 [ + - ][ + - ]: 2880 : stateProxy.ckLocalBranch()->insert( "DG", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
686 : : "comGhost" );
687 : :
688 : 5182 : auto d = Disc();
689 : 5182 : const auto& lid = d->Lid();
690 : 5182 : auto& inpofa = m_fd.Inpofa();
691 : 5182 : auto ncoord = m_coord[0].size();
692 : :
693 : : // nodelist with fromch, currently only used for an assert
694 : 5182 : [[maybe_unused]] const auto& nl = tk::cref_find( d->NodeCommMap(), fromch );
695 : :
696 : 5182 : auto& ghostelem = m_ghost[ fromch ]; // will associate to sender chare
697 : :
698 : : // Store ghost data coming from chare
699 [ + + ]: 111514 : for (const auto& g : ghost) { // loop over incoming ghost data
700 : 106332 : auto e = g.first; // remote/ghost tet id outside of chare boundary
701 : 106332 : const auto& nodes = std::get< 0 >( g.second ); // node IDs of face(s)
702 : 106332 : const auto& geo = std::get< 1 >( g.second ); // ghost elem geometry data
703 : 106332 : const auto& coordg = std::get< 2 >( g.second ); // coordinate of ghost node
704 : 106332 : const auto& inpoelg = std::get< 4 >( g.second ); // inpoel of ghost tet
705 : :
706 [ - + ][ - - ]: 106332 : Assert( nodes.size() % 3 == 0, "Face node IDs must be triplets" );
[ - - ][ - - ]
707 [ - + ][ - - ]: 106332 : Assert( nodes.size() <= 4*3, "Overflow of faces/tet received" );
[ - - ][ - - ]
708 [ - + ][ - - ]: 106332 : Assert( geo.size() % 5 == 0, "Ghost geometry size mismatch" );
[ - - ][ - - ]
709 [ - + ][ - - ]: 106332 : Assert( geo.size() == m_geoElem.nprop(), "Ghost geometry number mismatch" );
[ - - ][ - - ]
710 [ - + ][ - - ]: 106332 : Assert( coordg.size() == 3, "Incorrect ghost node coordinate size" );
[ - - ][ - - ]
711 [ - + ][ - - ]: 106332 : Assert( inpoelg.size() == 4, "Incorrect ghost inpoel size" );
[ - - ][ - - ]
712 : :
713 [ + + ]: 242562 : for (std::size_t n=0; n<nodes.size()/3; ++n) { // face(s) of ghost e
714 : : // node IDs of face on chare boundary
715 : 136230 : tk::UnsMesh::Face t{{ nodes[n*3+0], nodes[n*3+1], nodes[n*3+2] }};
716 : : // must find t in nodelist of chare-boundary adjacent to fromch
717 [ + - ][ + - ]: 136230 : Assert( nl.find(t[0]) != end(nl) &&
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
[ - - ]
718 : : nl.find(t[1]) != end(nl) &&
719 : : nl.find(t[2]) != end(nl),
720 : : "Ghost face not found in chare-node adjacency map on receiving end" );
721 : : // must find face in boundary-face adjacency map for fromch
722 [ + - ][ + - ]: 136230 : Assert( tk::cref_find(m_bndFace,fromch).find( t ) !=
[ + - ][ - + ]
[ - - ][ - - ]
[ - - ]
723 : : tk::cref_find(m_bndFace,fromch).cend(), "Ghost face not "
724 : : "found in boundary-face adjacency map on receiving end" );
725 : : // find local face & tet ids for t
726 [ + - ][ + - ]: 136230 : auto id = tk::cref_find( tk::cref_find(m_bndFace,fromch), t );
727 : : // compute face geometry for chare-boundary face
728 [ + - ]: 136230 : addGeoFace( t, id );
729 : : // add node-triplet to node-face connectivity
730 [ + - ]: 136230 : inpofa[3*id[0]+0] = tk::cref_find( lid, t[2] );
731 [ + - ]: 136230 : inpofa[3*id[0]+1] = tk::cref_find( lid, t[1] );
732 [ + - ]: 136230 : inpofa[3*id[0]+2] = tk::cref_find( lid, t[0] );
733 : :
734 : : // if ghost tet id not yet encountered on boundary with fromch
735 [ + - ]: 136230 : auto i = ghostelem.find( e );
736 [ + + ]: 136230 : if (i != end(ghostelem)) {
737 [ + - ]: 29898 : addEsuf( id, i->second ); // fill in elements surrounding face
738 [ + - ]: 29898 : addEsuel( id, i->second, t ); // fill in elements surrounding element
739 : : } else {
740 [ + - ]: 106332 : addEsuf( id, m_nunk ); // fill in elements surrounding face
741 [ + - ]: 106332 : addEsuel( id, m_nunk, t ); // fill in elements surrounding element
742 [ + - ]: 106332 : ghostelem[e] = m_nunk; // assign new local tet id to remote ghost id
743 [ + - ]: 106332 : m_geoElem.push_back( geo );// store ghost elem geometry
744 : 106332 : ++m_nunk; // increase number of unknowns on this chare
745 : 106332 : std::size_t counter = 0;
746 [ + + ]: 531660 : for (std::size_t gp=0; gp<4; ++gp) {
747 [ + - ]: 425328 : auto it = lid.find( inpoelg[gp] );
748 : : std::size_t lp;
749 [ + + ]: 425328 : if (it != end(lid))
750 : 357378 : lp = it->second;
751 : : else {
752 [ - + ][ - - ]: 67950 : Assert( nodes.size() == 3, "Expected node not found in lid" );
[ - - ][ - - ]
753 [ - + ][ - - ]: 67950 : Assert( gp == std::get< 3 >( g.second ),
[ - - ][ - - ]
754 : : "Ghost node not matching correct entry in ghost inpoel" );
755 : 67950 : lp = ncoord;
756 : 67950 : ++counter;
757 : : }
758 [ + - ]: 425328 : m_inpoel.push_back( lp ); // store ghost element connectivity
759 : : }
760 : : // only a single or no ghost node should be found
761 [ - + ][ - - ]: 106332 : Assert( counter <= 1, "Incorrect number of ghost nodes detected. "
[ - - ][ - - ]
[ - - ][ - - ]
762 : : "Detected "+ std::to_string(counter) +" ghost nodes" );
763 [ + + ]: 106332 : if (counter == 1) {
764 [ + - ]: 67950 : m_coord[0].push_back( coordg[0] ); // store ghost node coordinate
765 [ + - ]: 67950 : m_coord[1].push_back( coordg[1] );
766 [ + - ]: 67950 : m_coord[2].push_back( coordg[2] );
767 [ - + ][ - - ]: 67950 : Assert( m_inpoel[ 4*(m_nunk-1)+std::get< 3 >( g.second ) ] == ncoord,
[ - - ][ - - ]
768 : : "Mismatch in extended inpoel for ghost element" );
769 : 67950 : ++ncoord; // increase number of nodes on this chare
770 : : }
771 : : }
772 : :
773 : : // additional tests to ensure that entries in inpoel and t/inpofa match
774 [ + - ][ - + ]: 136230 : Assert( nodetripletMatch(id, t) == 3, "Mismatch/Overmatch in inpoel and "
[ - - ][ - - ]
[ - - ]
775 : : "inpofa at chare-boundary face" );
776 : : }
777 : : }
778 : :
779 : : // Signal the runtime system that all workers have received their
780 : : // face-adjacency
781 [ + + ]: 5182 : if (++m_nadj == m_ghostData.size()) faceAdj();
782 : 5182 : }
783 : :
784 : : std::size_t
785 : 136230 : DG::nodetripletMatch( const std::array< std::size_t, 2 >& id,
786 : : const tk::UnsMesh::Face& t )
787 : : // *****************************************************************************
788 : : // Check if entries in inpoel, inpofa and node-triplet are consistent
789 : : //! \param[in] id Local face and (inner) tet id adjacent to it
790 : : //! \param[in] t node-triplet associated with the chare boundary face
791 : : //! \return number of nodes in inpoel that matched with t and inpofa
792 : : // *****************************************************************************
793 : : {
794 : 136230 : const auto& lid = Disc()->Lid();
795 : 136230 : const auto& esuf = m_fd.Esuf();
796 : 136230 : const auto& inpofa = m_fd.Inpofa();
797 : :
798 : 136230 : std::size_t counter = 0;
799 [ + + ]: 681150 : for (std::size_t k=0; k<4; ++k)
800 : : {
801 : 544920 : auto el = esuf[ 2*id[0] ];
802 : 544920 : auto ip = m_inpoel[ 4*static_cast< std::size_t >( el )+k ];
803 [ - + ][ - - ]: 544920 : Assert( el == static_cast< int >( id[1] ), "Mismatch in id and esuf" );
[ - - ][ - - ]
804 [ + + ]: 2179680 : for (std::size_t j=0; j<3; ++j)
805 : : {
806 : 1634760 : auto jp = tk::cref_find( lid, t[j] );
807 : 1634760 : auto fp = inpofa[ 3*id[0]+(2-j) ];
808 [ + + ][ + - ]: 1634760 : if (ip == jp && ip == fp) ++counter;
809 : : }
810 : : }
811 : :
812 : 136230 : return counter;
813 : : }
814 : :
815 : : void
816 : 136230 : DG::addEsuf( const std::array< std::size_t, 2 >& id, std::size_t ghostid )
817 : : // *****************************************************************************
818 : : // Fill elements surrounding a face along chare boundary
819 : : //! \param[in] id Local face and (inner) tet id adjacent to it
820 : : //! \param[in] ghostid Local ID for ghost tet
821 : : //! \details This function extends and fills in the elements surrounding faces
822 : : //! data structure (esuf) so that the left and right element id is filled
823 : : //! in correctly on chare boundaries to contain the correct inner tet id and
824 : : //! the local tet id for the outer (ghost) tet, both adjacent to the given
825 : : //1 chare-face boundary. Prior to this function, this data structure does not
826 : : //! have yet face-element connectivity adjacent to chare-boundary faces, only
827 : : //! for physical boundaries and internal faces that are not on the chare
828 : : //! boundary (this latter purely as a result of mesh partitioning). The remote
829 : : //! element id of the ghost is stored in a location that is local to our own
830 : : //! esuf. The face numbering is such that esuf stores the element-face
831 : : //! connectivity first for the physical-boundary faces, followed by that of
832 : : //! the internal faces, followed by the chare-boundary faces. As a result,
833 : : //! esuf can be used by physics algorithms in exactly the same way as would be
834 : : //! used in serial. In serial, of course, this data structure is not extended
835 : : //! at the end by the chare-boundaries.
836 : : // *****************************************************************************
837 : : {
838 : 136230 : auto& esuf = m_fd.Esuf();
839 [ - + ][ - - ]: 136230 : Assert( 2*id[0]+1 < esuf.size(), "Indexing out of esuf" );
[ - - ][ - - ]
840 : :
841 : : // put in inner tet id
842 [ + - ][ + - ]: 136230 : Assert( esuf[ 2*id[0] ] == -2 && esuf[ 2*id[0]+1 ] == -2, "Updating esuf at "
[ - - ][ - - ]
[ - - ]
843 : : "wrong location instead of chare-boundary" );
844 : 136230 : esuf[ 2*id[0]+0 ] = static_cast< int >( id[1] );
845 : : // put in local id for outer/ghost tet
846 : 136230 : esuf[ 2*id[0]+1 ] = static_cast< int >( ghostid );
847 : 136230 : }
848 : :
849 : : void
850 : 136230 : DG::addEsuel( const std::array< std::size_t, 2 >& id,
851 : : std::size_t ghostid,
852 : : const tk::UnsMesh::Face& t )
853 : : // *****************************************************************************
854 : : // Fill elements surrounding a element along chare boundary
855 : : //! \param[in] id Local face and (inner) tet id adjacent to it
856 : : //! \param[in] ghostid Local ID for ghost tet
857 : : //! \param[in] t node-triplet associated with the chare boundary face
858 : : //! \details This function updates the elements surrounding element (esuel) data
859 : : // structure for the (inner) tets adjacent to the chare-boundaries. It fills
860 : : // esuel of this inner tet with the local tet-id that has been assigned to
861 : : // the outer ghost tet in DG::comGhost in place of the -1 before.
862 : : // *****************************************************************************
863 : : {
864 [ + - ]: 136230 : auto d = Disc();
865 : 136230 : [[maybe_unused]] const auto& esuf = m_fd.Esuf();
866 : 136230 : const auto& lid = d->Lid();
867 : :
868 : : std::array< tk::UnsMesh::Face, 4 > face;
869 [ + + ]: 681150 : for (std::size_t f = 0; f<4; ++f)
870 [ + + ]: 2179680 : for (std::size_t i = 0; i<3; ++i)
871 : 1634760 : face[f][i] = m_inpoel[ id[1]*4 + tk::lpofa[f][i] ];
872 : :
873 [ + - ]: 136230 : tk::UnsMesh::Face tl{{ tk::cref_find( lid, t[0] ),
874 [ + - ]: 136230 : tk::cref_find( lid, t[1] ),
875 [ + - ]: 272460 : tk::cref_find( lid, t[2] ) }};
876 : :
877 : 136230 : auto& esuel = m_fd.Esuel();
878 : 136230 : std::size_t i(0), nmatch(0);
879 [ + + ]: 681150 : for (const auto& f : face) {
880 [ + - ][ + + ]: 544920 : if (tk::UnsMesh::Eq< 3 >()( tl, f )) {
881 [ - + ][ - - ]: 136230 : Assert( esuel[ id[1]*4 + i ] == -1, "Incorrect boundary element found in "
[ - - ][ - - ]
882 : : "esuel");
883 : 136230 : esuel[ id[1]*4 + i ] = static_cast<int>(ghostid);
884 : 136230 : ++nmatch;
885 [ - + ][ - - ]: 136230 : Assert( esuel[ id[1]*4 + i ] == esuf[ 2*id[0]+1 ], "Incorrect boundary "
[ - - ][ - - ]
886 : : "element entered in esuel" );
887 [ - + ][ - - ]: 136230 : Assert( static_cast<int>(id[1]) == esuf[ 2*id[0]+0 ], "Boundary "
[ - - ][ - - ]
888 : : "element entered in incorrect esuel location" );
889 : : }
890 : 544920 : ++i;
891 : : }
892 : :
893 : : // ensure that exactly one face matched
894 [ - + ][ - - ]: 136230 : Assert( nmatch == 1, "Incorrect number of node-triplets (faces) matched for "
[ - - ][ - - ]
[ - - ]
895 : : "updating esuel; matching faces = "+ std::to_string(nmatch) );
896 : 136230 : }
897 : :
898 : : void
899 : 136230 : DG::addGeoFace( const tk::UnsMesh::Face& t,
900 : : const std::array< std::size_t, 2 >& id )
901 : : // *****************************************************************************
902 : : // Fill face-geometry data along chare boundary
903 : : //! \param[in] t Face (given by 3 global node IDs) on the chare boundary
904 : : //! \param[in] id Local face and (inner) tet id adjacent to face t
905 : : //! \details This function fills in the face geometry data along a chare
906 : : //! boundary.
907 : : // *****************************************************************************
908 : : {
909 [ + - ]: 136230 : auto d = Disc();
910 : 136230 : const auto& lid = d->Lid();
911 : :
912 : : // get global node IDs reversing order to get outward-pointing normal
913 [ + - ]: 136230 : auto A = tk::cref_find( lid, t[2] );
914 [ + - ]: 136230 : auto B = tk::cref_find( lid, t[1] );
915 [ + - ]: 136230 : auto C = tk::cref_find( lid, t[0] );
916 : 408690 : auto geochf = tk::geoFaceTri( {{m_coord[0][A], m_coord[0][B], m_coord[0][C]}},
917 : 408690 : {{m_coord[1][A], m_coord[1][B], m_coord[1][C]}},
918 [ + - ]: 1089840 : {{m_coord[2][A], m_coord[2][B], m_coord[2][C]}} );
919 : :
920 [ + + ]: 1089840 : for (std::size_t i=0; i<7; ++i)
921 [ + - ][ + - ]: 953610 : m_geoFace(id[0],i,0) = geochf(0,i,0);
922 : 136230 : }
923 : :
924 : : void
925 : 1026 : DG::faceAdj()
926 : : // *****************************************************************************
927 : : // Continue after face adjacency communication map completed on this chare
928 : : //! \details At this point the face communication map has been established
929 : : //! on this chare. Proceed to set up the nodal-comm map.
930 : : // *****************************************************************************
931 : : {
932 : 1026 : m_nadj = 0;
933 : :
934 : 1026 : tk::destroy(m_bndFace);
935 : :
936 : : // Ensure that all elements surrounding faces (are correct) including those at
937 : : // chare boundaries
938 [ + + ]: 1974186 : for (std::size_t f=0; f<m_nfac; ++f) {
939 [ - + ][ - - ]: 1973160 : Assert( m_fd.Esuf()[2*f] > -1,
[ - - ][ - - ]
940 : : "Left element in esuf cannot be physical ghost" );
941 [ + - ][ + + ]: 1973160 : if (f >= m_fd.Nbfac())
942 [ - + ][ - - ]: 1729994 : Assert( m_fd.Esuf()[2*f+1] > -1,
[ - - ][ - - ]
943 : : "Right element in esuf for internal/chare faces cannot be a ghost" );
944 : : }
945 : :
946 : : // Ensure that all elements surrounding elements are correct including those
947 : : // at chare boundaries
948 : 1026 : const auto& esuel = m_fd.Esuel();
949 : 1026 : std::size_t nbound = 0;
950 [ + + ]: 892757 : for (std::size_t e=0; e<esuel.size()/4; ++e) {
951 [ + + ]: 4458655 : for (std::size_t f=0; f<4; ++f)
952 [ + + ]: 3566924 : if (esuel[4*e+f] == -1) ++nbound;
953 : : }
954 [ + - ][ - + ]: 1026 : Assert( nbound == m_fd.Nbfac(), "Incorrect number of ghost-element -1's in "
[ - - ][ - - ]
[ - - ]
955 : : "updated esuel" );
956 : :
957 : : // Error checking on ghost data
958 [ + + ]: 6208 : for(const auto& n : m_ghostData)
959 [ + + ]: 111514 : for([[maybe_unused]] const auto& i : n.second)
960 [ - + ][ - - ]: 106332 : Assert( i.first < m_fd.Esuel().size()/4, "Sender contains ghost tet id " );
[ - - ][ - - ]
961 : :
962 : : // Perform leak test on face geometry data structure enlarged by ghosts
963 [ + - ][ - + ]: 1026 : Assert( !leakyAdjacency(), "Face adjacency leaky" );
[ - - ][ - - ]
[ - - ]
964 [ + - ][ - + ]: 1026 : Assert( faceMatch(), "Chare-boundary element-face connectivity (esuf) does "
[ - - ][ - - ]
[ - - ]
965 : : "not match" );
966 : :
967 : : // Create new map of elements along chare boundary which are ghosts for
968 : : // neighboring chare, associated with that chare ID
969 [ + + ]: 6208 : for (const auto& [cid, cgd] : m_ghostData)
970 : : {
971 [ + - ]: 5182 : auto& sg = m_sendGhost[cid];
972 [ + + ]: 111514 : for (const auto& e : cgd)
973 : : {
974 [ + - ][ - + ]: 106332 : Assert(sg.find(e.first) == sg.end(), "Repeating element found in "
[ - - ][ - - ]
[ - - ]
975 : : "ghost data");
976 [ + - ]: 106332 : sg.insert(e.first);
977 : : }
978 [ - + ][ - - ]: 5182 : Assert(sg.size() == cgd.size(), "Incorrect size for sendGhost");
[ - - ][ - - ]
979 : : }
980 [ - + ][ - - ]: 1026 : Assert(m_sendGhost.size() == m_ghostData.size(), "Incorrect number of "
[ - - ][ - - ]
981 : : "chares in sendGhost");
982 : :
983 : : // Error checking on ghost data
984 [ + + ]: 6208 : for(const auto& n : m_sendGhost)
985 [ + + ]: 111514 : for([[maybe_unused]] const auto& i : n.second)
986 [ - + ][ - - ]: 106332 : Assert( i < m_fd.Esuel().size()/4, "Sender contains ghost tet id. " );
[ - - ][ - - ]
987 : :
988 : : // Generate and store Esup data-structure in a map
989 [ + - ]: 1026 : auto esup = tk::genEsup(m_inpoel, 4);
990 [ + - ][ + + ]: 254759 : for (std::size_t p=0; p<Disc()->Gid().size(); ++p)
991 : : {
992 [ + + ]: 4178035 : for (auto e : tk::Around(esup, p))
993 : : {
994 : : // since inpoel has been augmented with the face-ghost cell previously,
995 : : // esup also contains cells which are not on this mesh-chunk, hence the
996 : : // following test
997 [ + + ][ + - ]: 3924302 : if (e < m_fd.Esuel().size()/4) m_esup[p].push_back(e);
[ + - ]
998 : : }
999 : : }
1000 : :
1001 : : // Error checking on Esup map
1002 [ + + ]: 254759 : for(const auto& p : m_esup)
1003 [ + + ]: 3820657 : for([[maybe_unused]] const auto& e : p.second)
1004 [ - + ][ - - ]: 3566924 : Assert( e < m_fd.Esuel().size()/4, "Esup contains tet id greater than "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
1005 : : + std::to_string(m_fd.Esuel().size()/4-1) +" : "+ std::to_string(e) );
1006 : :
1007 [ + - ]: 1026 : auto meshid = Disc()->MeshId();
1008 [ + - ]: 1026 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1009 [ + - ][ + - ]: 2052 : CkCallback(CkReductionTarget(Transporter,startEsup), Disc()->Tr()) );
[ + - ]
1010 : 1026 : }
1011 : :
1012 : : void
1013 : 1026 : DG::nodeNeighSetup()
1014 : : // *****************************************************************************
1015 : : // Setup node-neighborhood (esup)
1016 : : //! \details At this point the face-ghost communication map has been established
1017 : : //! on this chare. This function begins generating the node-ghost comm map.
1018 : : // *****************************************************************************
1019 : : {
1020 [ + + ]: 1026 : if (Disc()->NodeCommMap().empty())
1021 : : // in serial, skip setting up node-neighborhood
1022 : 43 : { comesup_complete(); }
1023 : : else
1024 : : {
1025 : 983 : const auto& nodeCommMap = Disc()->NodeCommMap();
1026 : :
1027 : : // send out node-neighborhood map
1028 [ + + ]: 9211 : for (const auto& [cid, nlist] : nodeCommMap)
1029 : : {
1030 : 16456 : std::unordered_map< std::size_t, std::vector< std::size_t > > bndEsup;
1031 : 8228 : std::unordered_map< std::size_t, std::vector< tk::real > > nodeBndCells;
1032 [ + + ]: 133890 : for (const auto& p : nlist)
1033 : : {
1034 [ + - ][ + - ]: 125662 : auto pl = tk::cref_find(Disc()->Lid(), p);
1035 : : // fill in the esup for the chare-boundary
1036 [ + - ]: 125662 : const auto& pesup = tk::cref_find(m_esup, pl);
1037 [ + - ][ + - ]: 125662 : bndEsup[p] = pesup;
1038 : :
1039 : : // fill a map with the element ids from esup as keys and geoElem as
1040 : : // values, and another map containing these elements associated with
1041 : : // the chare id with which they are node-neighbors.
1042 [ + + ]: 1107729 : for (const auto& e : pesup)
1043 : : {
1044 [ + - ][ + - ]: 982067 : nodeBndCells[e] = m_geoElem[e];
1045 : :
1046 : : // add these esup-elements into map of elements along chare boundary
1047 [ - + ][ - - ]: 982067 : Assert( e < m_fd.Esuel().size()/4, "Sender contains ghost tet id." );
[ - - ][ - - ]
1048 [ + - ][ + - ]: 982067 : m_sendGhost[cid].insert(e);
1049 : : }
1050 : : }
1051 : :
1052 [ + - ][ + - ]: 8228 : thisProxy[cid].comEsup(thisIndex, bndEsup, nodeBndCells);
1053 : : }
1054 : : }
1055 : :
1056 : 1026 : ownesup_complete();
1057 : 1026 : }
1058 : :
1059 : : void
1060 : 8228 : DG::comEsup( int fromch,
1061 : : const std::unordered_map< std::size_t, std::vector< std::size_t > >& bndEsup,
1062 : : const std::unordered_map< std::size_t, std::vector< tk::real > >&
1063 : : nodeBndCells )
1064 : : // *****************************************************************************
1065 : : //! \brief Receive elements-surrounding-points data-structure for points on
1066 : : // common boundary between receiving and sending neighbor chare, and the
1067 : : // element geometries for these new elements
1068 : : //! \param[in] fromch Sender chare id
1069 : : //! \param[in] bndEsup Elements-surrounding-points data-structure from fromch
1070 : : //! \param[in] nodeBndCells Map containing element geometries associated with
1071 : : //! remote element IDs in the esup
1072 : : // *****************************************************************************
1073 : : {
1074 : 8228 : auto& chghost = m_ghost[fromch];
1075 : :
1076 : : // Extend remote-local element id map and element geometry array
1077 [ + + ]: 488199 : for (const auto& e : nodeBndCells)
1078 : : {
1079 : : // need to check following, because 'e' could have been added previously in
1080 : : // remote-local element id map as a part of face-communication, i.e. as a
1081 : : // face-ghost element
1082 [ + - ][ + + ]: 479971 : if (chghost.find(e.first) == chghost.end())
1083 : : {
1084 [ + - ]: 373639 : chghost[e.first] = m_nunk;
1085 [ + - ]: 373639 : m_geoElem.push_back(e.second);
1086 : 373639 : ++m_nunk;
1087 : : }
1088 : : }
1089 : :
1090 : : // Store incoming data in comm-map buffer for Esup
1091 [ + + ]: 133890 : for (const auto& [node, elist] : bndEsup)
1092 : : {
1093 [ + - ][ + - ]: 125662 : auto pl = tk::cref_find(Disc()->Lid(), node);
1094 [ + - ]: 125662 : auto& pesup = m_esupc[pl];
1095 [ + + ]: 1107729 : for (auto e : elist)
1096 : : {
1097 [ + - ]: 982067 : auto el = tk::cref_find(chghost, e);
1098 [ + - ]: 982067 : pesup.push_back(el);
1099 : : }
1100 : : }
1101 : :
1102 : : // if we have heard from all fellow chares that we share at least a single
1103 : : // node, edge, or face with
1104 [ + + ]: 8228 : if (++m_ncomEsup == Disc()->NodeCommMap().size()) {
1105 : 983 : m_ncomEsup = 0;
1106 : 983 : comesup_complete();
1107 : : }
1108 : 8228 : }
1109 : :
1110 : : void
1111 : 1026 : DG::adj()
1112 : : // *****************************************************************************
1113 : : // Finish up with adjacency maps, and do a global-sync to begin problem setup
1114 : : //! \details At this point, the nodal- and face-adjacency has been set up. This
1115 : : // function does some error checking on the nodal-adjacency and prepares
1116 : : // for problem setup.
1117 : : // *****************************************************************************
1118 : : {
1119 : : // combine own and communicated contributions to elements surrounding points
1120 [ + + ]: 85882 : for (auto& [p, elist] : m_esupc)
1121 : : {
1122 [ + - ]: 84856 : auto& pesup = tk::ref_find(m_esup, p);
1123 [ + + ]: 1066923 : for ([[maybe_unused]] auto e : elist)
1124 : : {
1125 [ - + ][ - - ]: 982067 : Assert( e >= m_fd.Esuel().size()/4, "Non-ghost element received from "
[ - - ][ - - ]
1126 : : "esup buffer." );
1127 : : }
1128 [ + - ]: 84856 : tk::concat< std::size_t >(std::move(elist), pesup);
1129 : : }
1130 : :
1131 : 1026 : tk::destroy(m_ghostData);
1132 : 1026 : tk::destroy(m_esupc);
1133 : :
1134 [ - + ][ - - ]: 1026 : if ( g_inputdeck.get< tag::cmd, tag::feedback >() ) Disc()->Tr().chadj();
[ - - ]
1135 : :
1136 : : // Error checking on ghost data
1137 [ + + ]: 9254 : for(const auto& n : m_sendGhost)
1138 [ + + ]: 488199 : for([[maybe_unused]] const auto& i : n.second)
1139 [ - + ][ - - ]: 479971 : Assert( i < m_fd.Esuel().size()/4, "Sender contains ghost tet id. ");
[ - - ][ - - ]
1140 : :
1141 : : // Resize solution vectors, lhs and rhs by the number of ghost tets
1142 [ + - ]: 1026 : m_u.resize( m_nunk );
1143 [ + - ]: 1026 : m_un.resize( m_nunk );
1144 [ + - ]: 1026 : m_p.resize( m_nunk );
1145 [ + - ]: 1026 : m_lhs.resize( m_nunk );
1146 [ + - ]: 1026 : m_rhs.resize( m_nunk );
1147 : :
1148 : : // Create a mapping between local ghost tet ids and zero-based boundary ids
1149 [ + - ]: 2052 : std::vector< std::size_t > c( tk::sumvalsize( m_ghost ) );
1150 : 1026 : std::size_t j = 0;
1151 [ + + ]: 9254 : for (const auto& n : m_ghost) {
1152 [ + + ]: 488199 : for(const auto& i : n.second) {
1153 : 479971 : c[j++] = i.second;
1154 : : }
1155 : : }
1156 [ + - ]: 1026 : m_bid = tk::assignLid( c );
1157 : :
1158 : : // Size communication buffer that receives number of degrees of freedom
1159 [ + + ][ + - ]: 4104 : for (auto& n : m_ndofc) n.resize( m_bid.size() );
1160 [ + + ][ + - ]: 4104 : for (auto& u : m_uc) u.resize( m_bid.size() );
1161 [ + + ][ + - ]: 4104 : for (auto& p : m_pc) p.resize( m_bid.size() );
1162 : :
1163 : : // Initialize number of degrees of freedom in mesh elements
1164 : 1026 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
1165 [ + + ]: 1026 : if( pref )
1166 : : {
1167 : 134 : const auto ndofmax = g_inputdeck.get< tag::pref, tag::ndofmax >();
1168 [ + - ]: 134 : m_ndof.resize( m_nunk, ndofmax );
1169 : : }
1170 : : else
1171 : : {
1172 : 892 : const auto ndof = g_inputdeck.get< tag::discr, tag::ndof >();
1173 [ + - ]: 892 : m_ndof.resize( m_nunk, ndof );
1174 : : }
1175 : :
1176 : : // Ensure that we also have all the geometry and connectivity data
1177 : : // (including those of ghosts)
1178 [ - + ][ - - ]: 1026 : Assert( m_geoElem.nunk() == m_u.nunk(), "GeoElem unknowns size mismatch" );
[ - - ][ - - ]
1179 : :
1180 : : // Basic error checking on ghost tet ID map
1181 [ + - ][ - + ]: 1026 : Assert( m_ghost.find( thisIndex ) == m_ghost.cend(),
[ - - ][ - - ]
[ - - ]
1182 : : "Ghost id map should not contain data for own chare ID" );
1183 : :
1184 : : // Store expected ghost tet IDs
1185 [ + + ]: 9254 : for (const auto& n : m_ghost)
1186 [ + + ]: 488199 : for ([[maybe_unused]] const auto& g : n.second)
1187 [ + - ][ - + ]: 479971 : Assert( m_exptGhost.insert( g.second ).second,
[ - - ][ - - ]
[ - - ]
1188 : : "Failed to store local tetid as exptected ghost id" );
1189 : :
1190 : : // Signal the runtime system that all workers have received their adjacency
1191 [ + - ][ + - ]: 1026 : std::vector< std::size_t > meshdata{ m_initial, Disc()->MeshId() };
1192 [ + - ]: 1026 : contribute( meshdata, CkReduction::sum_ulong,
1193 [ + - ][ + - ]: 2052 : CkCallback(CkReductionTarget(Transporter,comfinal), Disc()->Tr()) );
[ + - ]
1194 : 1026 : }
1195 : :
1196 : : void
1197 : 666 : DG::registerReducers()
1198 : : // *****************************************************************************
1199 : : // Configure Charm++ reduction types
1200 : : //! \details Since this is a [initnode] routine, the runtime system executes the
1201 : : //! routine exactly once on every logical node early on in the Charm++ init
1202 : : //! sequence. Must be static as it is called without an object. See also:
1203 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
1204 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
1205 : : // *****************************************************************************
1206 : : {
1207 : 666 : ElemDiagnostics::registerReducers();
1208 : 666 : }
1209 : :
1210 : : void
1211 : 20822 : DG::ResumeFromSync()
1212 : : // *****************************************************************************
1213 : : // Return from migration
1214 : : //! \details This is called when load balancing (LB) completes. The presence of
1215 : : //! this function does not affect whether or not we block on LB.
1216 : : // *****************************************************************************
1217 : : {
1218 [ - + ][ - - ]: 20822 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
1219 : :
1220 [ + - ]: 20822 : if (!g_inputdeck.get< tag::cmd, tag::nonblocking >()) next();
1221 : 20822 : }
1222 : :
1223 : : void
1224 : 918 : DG::setup()
1225 : : // *****************************************************************************
1226 : : // Set initial conditions, generate lhs, output mesh
1227 : : // *****************************************************************************
1228 : : {
1229 [ + + ][ + + ]: 1798 : if (g_inputdeck.get< tag::cmd, tag::chare >() ||
1230 [ + + ]: 880 : g_inputdeck.get< tag::cmd, tag::quiescence >())
1231 [ + - ][ + - ]: 541 : stateProxy.ckLocalBranch()->insert( "DG", thisIndex, CkMyPe(), Disc()->It(),
[ + - ][ + - ]
1232 : : "setup" );
1233 : :
1234 : 918 : auto d = Disc();
1235 : :
1236 : : // Basic error checking on sizes of element geometry data and connectivity
1237 [ - + ][ - - ]: 918 : Assert( m_geoElem.nunk() == m_lhs.nunk(), "Size mismatch in DG::setup()" );
[ - - ][ - - ]
1238 : :
1239 : : // Compute left-hand side of discrete PDEs
1240 : 918 : lhs();
1241 : :
1242 : : // Determine elements inside user-defined IC box
1243 [ + + ]: 1836 : for (auto& eq : g_dgpde)
1244 [ + - ]: 918 : eq.IcBoxElems( m_geoElem, m_fd.Esuel().size()/4, m_boxelems );
1245 : :
1246 : : // Compute volume of user-defined box IC
1247 [ + - ]: 918 : d->boxvol( {} ); // punt for now
1248 : :
1249 : : // Query time history field output labels from all PDEs integrated
1250 : 918 : const auto& hist_points = g_inputdeck.get< tag::history, tag::point >();
1251 [ - + ]: 918 : if (!hist_points.empty()) {
1252 : 0 : std::vector< std::string > histnames;
1253 [ - - ]: 0 : for (const auto& eq : g_dgpde) {
1254 [ - - ]: 0 : auto n = eq.histNames();
1255 [ - - ]: 0 : histnames.insert( end(histnames), begin(n), end(n) );
1256 : : }
1257 [ - - ]: 0 : d->histheader( std::move(histnames) );
1258 : : }
1259 : 918 : }
1260 : :
1261 : : void
1262 : 918 : DG::box( tk::real v )
1263 : : // *****************************************************************************
1264 : : // Receive total box IC volume and set conditions in box
1265 : : //! \param[in] v Total volume within user-specified box
1266 : : // *****************************************************************************
1267 : : {
1268 : 918 : auto d = Disc();
1269 : :
1270 : : // Store user-defined box IC volume
1271 : 918 : d->Boxvol() = v;
1272 : :
1273 : : // Set initial conditions for all PDEs
1274 [ + + ]: 1836 : for (const auto& eq : g_dgpde)
1275 : : {
1276 [ + - ]: 918 : eq.initialize( m_lhs, m_inpoel, m_coord, m_boxelems, m_u, d->T(),
1277 : 918 : m_fd.Esuel().size()/4 );
1278 [ + - ]: 918 : eq.updatePrimitives( m_u, m_lhs, m_geoElem, m_p, m_fd.Esuel().size()/4 );
1279 : : }
1280 : :
1281 : 918 : m_un = m_u;
1282 : :
1283 : : // Output initial conditions to file (regardless of whether it was requested)
1284 [ + - ][ + - ]: 918 : startFieldOutput( CkCallback(CkIndex_DG::start(), thisProxy[thisIndex]) );
[ + - ]
1285 : 918 : }
1286 : :
1287 : : void
1288 : 918 : DG::start()
1289 : : // *****************************************************************************
1290 : : // Start time stepping
1291 : : // *****************************************************************************
1292 : : {
1293 : : // Free memory storing output mesh
1294 : 918 : m_outmesh.destroy();
1295 : :
1296 : : // Start timer measuring time stepping wall clock time
1297 : 918 : Disc()->Timer().zero();
1298 : : // Zero grind-timer
1299 : 918 : Disc()->grindZero();
1300 : : // Start time stepping by computing the size of the next time step)
1301 : 918 : next();
1302 : 918 : }
1303 : :
1304 : : void
1305 : 25778 : DG::startFieldOutput( CkCallback c )
1306 : : // *****************************************************************************
1307 : : // Start preparing fields for output to file
1308 : : //! \param[in] c Function to continue with after the write
1309 : : // *****************************************************************************
1310 : : {
1311 : : // No field output in benchmark mode or if field output frequency not hit
1312 [ + + ][ + + ]: 25778 : if (g_inputdeck.get< tag::cmd, tag::benchmark >() || !fieldOutput()) {
[ + + ]
1313 : :
1314 : 23014 : c.send();
1315 : :
1316 : : } else {
1317 : :
1318 : : // Optionally refine mesh for field output
1319 : 2764 : auto d = Disc();
1320 : :
1321 [ + + ]: 2764 : if (refinedOutput()) {
1322 : :
1323 [ + - ]: 33 : const auto& tr = tk::remap( m_fd.Triinpoel(), d->Gid() );
1324 [ + - ][ + - ]: 33 : d->Ref()->outref( m_fd.Bface(), {}, tr, c );
1325 : :
1326 : : } else {
1327 : :
1328 : : // cut off ghosts from mesh connectivity and coordinates
1329 [ + - ]: 2731 : const auto& tr = tk::remap( m_fd.Triinpoel(), d->Gid() );
1330 [ + - ]: 2731 : extractFieldOutput( {}, d->Chunk(), d->Coord(), {}, {},
1331 : : d->NodeCommMap(), m_fd.Bface(), {}, tr, c );
1332 : :
1333 : : }
1334 : :
1335 : : }
1336 : 25778 : }
1337 : :
1338 : : void
1339 : 74580 : DG::next()
1340 : : // *****************************************************************************
1341 : : // Advance equations to next time step
1342 : : // *****************************************************************************
1343 : : {
1344 : 74580 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
1345 : :
1346 : 74580 : auto d = Disc();
1347 : :
1348 [ + + ][ + + ]: 74580 : if (pref && m_stage == 0 && d->T() > 0)
[ + + ][ + + ]
1349 [ + + ]: 3272 : for (const auto& eq : g_dgpde)
1350 [ + - ]: 1636 : eq.eval_ndof( m_nunk, m_coord, m_inpoel, m_fd, m_u,
1351 : 1636 : g_inputdeck.get< tag::pref, tag::indicator >(),
1352 : 1636 : g_inputdeck.get< tag::discr, tag::ndof >(),
1353 : 1636 : g_inputdeck.get< tag::pref, tag::ndofmax >(),
1354 : 1636 : g_inputdeck.get< tag::pref, tag::tolref >(),
1355 : 1636 : m_ndof );
1356 : :
1357 : : // communicate solution ghost data (if any)
1358 [ + + ]: 74580 : if (m_sendGhost.empty())
1359 : 4395 : comsol_complete();
1360 : : else
1361 [ + + ]: 569715 : for(const auto& [cid, ghostdata] : m_sendGhost) {
1362 [ + - ]: 999060 : std::vector< std::size_t > tetid( ghostdata.size() );
1363 [ + - ]: 999060 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
1364 [ + - ]: 999060 : prim( ghostdata.size() );
1365 : 499530 : std::vector< std::size_t > ndof;
1366 : 499530 : std::size_t j = 0;
1367 [ + + ]: 19317180 : for(const auto& i : ghostdata) {
1368 [ - + ][ - - ]: 18817650 : Assert( i < m_fd.Esuel().size()/4, "Sending solution ghost data" );
[ - - ][ - - ]
1369 : 18817650 : tetid[j] = i;
1370 [ + - ]: 18817650 : u[j] = m_u[i];
1371 [ + - ]: 18817650 : prim[j] = m_p[i];
1372 [ + + ][ + + ]: 18817650 : if (pref && m_stage == 0) ndof.push_back( m_ndof[i] );
[ + - ]
1373 : 18817650 : ++j;
1374 : : }
1375 [ + - ][ + - ]: 499530 : thisProxy[ cid ].comsol( thisIndex, m_stage, tetid, u, prim, ndof );
1376 : : }
1377 : :
1378 : 74580 : ownsol_complete();
1379 : 74580 : }
1380 : :
1381 : : void
1382 : 499530 : DG::comsol( int fromch,
1383 : : std::size_t fromstage,
1384 : : const std::vector< std::size_t >& tetid,
1385 : : const std::vector< std::vector< tk::real > >& u,
1386 : : const std::vector< std::vector< tk::real > >& prim,
1387 : : const std::vector< std::size_t >& ndof )
1388 : : // *****************************************************************************
1389 : : // Receive chare-boundary solution ghost data from neighboring chares
1390 : : //! \param[in] fromch Sender chare id
1391 : : //! \param[in] fromstage Sender chare time step stage
1392 : : //! \param[in] tetid Ghost tet ids we receive solution data for
1393 : : //! \param[in] u Solution ghost data
1394 : : //! \param[in] prim Primitive variables in ghost cells
1395 : : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
1396 : : //! \details This function receives contributions to the unlimited solution
1397 : : //! from fellow chares.
1398 : : // *****************************************************************************
1399 : : {
1400 [ - + ][ - - ]: 499530 : Assert( u.size() == tetid.size(), "Size mismatch in DG::comsol()" );
[ - - ][ - - ]
1401 [ - + ][ - - ]: 499530 : Assert( prim.size() == tetid.size(), "Size mismatch in DG::comsol()" );
[ - - ][ - - ]
1402 : :
1403 : 499530 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
1404 : :
1405 [ + + ][ + + ]: 499530 : if (pref && fromstage == 0)
1406 [ - + ][ - - ]: 13280 : Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comsol()" );
[ - - ][ - - ]
1407 : :
1408 : : // Find local-to-ghost tet id map for sender chare
1409 : 499530 : const auto& n = tk::cref_find( m_ghost, fromch );
1410 : :
1411 [ + + ]: 19317180 : for (std::size_t i=0; i<tetid.size(); ++i) {
1412 [ + - ]: 18817650 : auto j = tk::cref_find( n, tetid[i] );
1413 [ - + ][ - - ]: 18817650 : Assert( j >= m_fd.Esuel().size()/4, "Receiving solution non-ghost data" );
[ - - ][ - - ]
1414 [ + - ]: 18817650 : auto b = tk::cref_find( m_bid, j );
1415 [ - + ][ - - ]: 18817650 : Assert( b < m_uc[0].size(), "Indexing out of bounds" );
[ - - ][ - - ]
1416 [ + - ]: 18817650 : m_uc[0][b] = u[i];
1417 [ + - ]: 18817650 : m_pc[0][b] = prim[i];
1418 [ + + ][ + + ]: 18817650 : if (pref && fromstage == 0) {
1419 [ - + ][ - - ]: 395110 : Assert( b < m_ndofc[0].size(), "Indexing out of bounds" );
[ - - ][ - - ]
1420 : 395110 : m_ndofc[0][b] = ndof[i];
1421 : : }
1422 : : }
1423 : :
1424 : : // if we have received all solution ghost contributions from neighboring
1425 : : // chares (chares we communicate along chare-boundary faces with), and
1426 : : // contributed our solution to these neighbors, proceed to reconstructions
1427 [ + + ]: 499530 : if (++m_nsol == m_sendGhost.size()) {
1428 : 70185 : m_nsol = 0;
1429 : 70185 : comsol_complete();
1430 : : }
1431 : 499530 : }
1432 : :
1433 : : void
1434 : 2764 : DG::extractFieldOutput(
1435 : : const std::vector< std::size_t >& /*ginpoel*/,
1436 : : const tk::UnsMesh::Chunk& chunk,
1437 : : const tk::UnsMesh::Coords& coord,
1438 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /*addedNodes*/,
1439 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
1440 : : const tk::NodeCommMap& nodeCommMap,
1441 : : const std::map< int, std::vector< std::size_t > >& bface,
1442 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
1443 : : const std::vector< std::size_t >& triinpoel,
1444 : : CkCallback c )
1445 : : // *****************************************************************************
1446 : : // Extract field output going to file
1447 : : //! \param[in] chunk Field-output mesh chunk (connectivity and global<->local
1448 : : //! id maps)
1449 : : //! \param[in] coord Field-output mesh node coordinates
1450 : : //! \param[in] addedTets Field-output mesh cells and their parents (local ids)
1451 : : //! \param[in] nodeCommMap Field-output mesh node communication map
1452 : : //! \param[in] bface Field-output meshndary-faces mapped to side set ids
1453 : : //! \param[in] triinpoel Field-output mesh boundary-face connectivity
1454 : : //! \param[in] c Function to continue with after the write
1455 : : // *****************************************************************************
1456 : : {
1457 [ + - ]: 2764 : m_outmesh.chunk = chunk;
1458 [ + - ]: 2764 : m_outmesh.coord = coord;
1459 [ + - ]: 2764 : m_outmesh.triinpoel = triinpoel;
1460 [ + - ]: 2764 : m_outmesh.bface = bface;
1461 [ + - ]: 2764 : m_outmesh.nodeCommMap = nodeCommMap;
1462 : :
1463 : 2764 : const auto& inpoel = std::get< 0 >( chunk );
1464 : 2764 : auto nelem = inpoel.size() / 4;
1465 : :
1466 : : // Evaluate element solution on incoming mesh
1467 [ + - ]: 5528 : auto [ue,pe,un,pn] = evalSolution( inpoel, coord, addedTets );
1468 : :
1469 : : // Collect field output from numerical solution requested by user
1470 [ + - ]: 2764 : m_elemfields = numericFieldOutput( ue, tk::Centering::ELEM, pe );
1471 [ + - ]: 2764 : m_nodefields = numericFieldOutput( un, tk::Centering::NODE, pn );
1472 : :
1473 : : // Collect field output from analytical solutions (if exist)
1474 [ + - ]: 5528 : auto geoElem = tk::genGeoElemTet( inpoel, coord );
1475 [ + - ]: 2764 : auto t = Disc()->T();
1476 [ + + ]: 5528 : for (const auto& eq : g_dgpde) {
1477 [ + - ][ + - ]: 2764 : analyticFieldOutput( eq, tk::Centering::ELEM, geoElem.extract(1,0),
1478 [ + - ][ + - ]: 5528 : geoElem.extract(2,0), geoElem.extract(3,0), t, m_elemfields );
1479 [ + - ]: 2764 : analyticFieldOutput( eq, tk::Centering::NODE, coord[0], coord[1], coord[2],
1480 : 2764 : t, m_nodefields );
1481 : : }
1482 : :
1483 : : // Add adaptive indicator array to element-centered field output
1484 [ + + ]: 2764 : if (g_inputdeck.get< tag::pref, tag::pref >()) {
1485 [ + - ]: 804 : std::vector< tk::real > ndof( begin(m_ndof), end(m_ndof) );
1486 [ + - ]: 402 : ndof.resize( nelem );
1487 [ + + ]: 87834 : for (const auto& [child,parent] : addedTets)
1488 : 87432 : ndof[child] = static_cast< tk::real >( m_ndof[parent] );
1489 [ + - ]: 402 : m_elemfields.push_back( ndof );
1490 : : }
1491 : :
1492 : : // Add shock detection marker array to element-centered field output
1493 [ + - ]: 2764 : std::vector< tk::real > shockmarker( begin(m_shockmarker), end(m_shockmarker) );
1494 : : // Here m_shockmarker has a size of m_u.nunk() which is the number of the
1495 : : // elements within this partition (nelem) plus the ghost partition cells. In
1496 : : // terms of output purpose, we only need the solution data within this
1497 : : // partition. Therefore, resizing it to nelem removes the extra partition
1498 : : // boundary allocations in the shockmarker vector. Since the code assumes that
1499 : : // the boundary elements are on the top, the resize operation keeps the lower
1500 : : // portion.
1501 [ + - ]: 2764 : shockmarker.resize( nelem );
1502 [ + + ]: 160276 : for (const auto& [child,parent] : addedTets)
1503 : 157512 : shockmarker[child] = static_cast< tk::real >(m_shockmarker[parent]);
1504 [ + - ]: 2764 : m_elemfields.push_back( shockmarker );
1505 : :
1506 : : // Send node fields contributions to neighbor chares
1507 [ + + ]: 2764 : if (nodeCommMap.empty())
1508 [ + - ]: 200 : comnodeout_complete();
1509 : : else {
1510 : 2564 : const auto& lid = std::get< 2 >( chunk );
1511 [ + - ]: 5128 : auto esup = tk::genEsup( inpoel, 4 );
1512 [ + + ]: 21250 : for(const auto& [ch,nodes] : nodeCommMap) {
1513 : : // Pack node field data in chare boundary nodes
1514 : : std::vector< std::vector< tk::real > >
1515 [ + - ][ + - ]: 56058 : l( m_nodefields.size(), std::vector< tk::real >( nodes.size() ) );
1516 [ - + ]: 18686 : for (std::size_t f=0; f<m_nodefields.size(); ++f) {
1517 : 0 : std::size_t j = 0;
1518 [ - - ]: 0 : for (auto g : nodes)
1519 [ - - ]: 0 : l[f][j++] = m_nodefields[f][ tk::cref_find(lid,g) ];
1520 : : }
1521 : : // Pack (partial) number of elements surrounding chare boundary nodes
1522 [ + - ]: 18686 : std::vector< std::size_t > nesup( nodes.size() );
1523 : 18686 : std::size_t j = 0;
1524 [ + + ]: 179040 : for (auto g : nodes) {
1525 [ + - ]: 160354 : auto i = tk::cref_find( lid, g );
1526 : 160354 : nesup[j++] = esup.second[i+1] - esup.second[i];
1527 : : }
1528 [ + - ][ + - ]: 56058 : thisProxy[ch].comnodeout(
1529 [ + - ]: 37372 : std::vector<std::size_t>(begin(nodes),end(nodes)), nesup, l );
1530 : : }
1531 : : }
1532 : :
1533 [ + - ]: 2764 : ownnod_complete( c );
1534 : 2764 : }
1535 : :
1536 : : std::tuple< tk::Fields, tk::Fields, tk::Fields, tk::Fields >
1537 : 2764 : DG::evalSolution(
1538 : : const std::vector< std::size_t >& inpoel,
1539 : : const tk::UnsMesh::Coords& coord,
1540 : : const std::unordered_map< std::size_t, std::size_t >& addedTets )
1541 : : // *****************************************************************************
1542 : : // Evaluate solution on incomping (a potentially refined) mesh
1543 : : //! \param[in] inpoel Incoming (potentially refined field-output) mesh
1544 : : //! connectivity
1545 : : //! \param[in] coord Incoming (potentially refined Field-output) mesh node
1546 : : //! coordinates
1547 : : //! \param[in] addedTets Field-output mesh cells and their parents (local ids)
1548 : : //! \details This function evaluates the solution on the incoming mesh. The
1549 : : //! incoming mesh can be refined but can also be just the mesh the numerical
1550 : : //! solution is computed on.
1551 : : //! \note If the incoming mesh is refined (for field putput) compared to the
1552 : : //! mesh the numerical solution is computed on, the solution is evaluated in
1553 : : //! cells as wells as in nodes. If the solution is not refined, the solution
1554 : : //! is evaluated in nodes.
1555 : : //! \return Solution in cells, primitive variables in cells, solution in nodes,
1556 : : //! primitive variables in nodes of incoming mesh.
1557 : : // *****************************************************************************
1558 : : {
1559 : : using tk::dot;
1560 : : using tk::real;
1561 : :
1562 : 2764 : const auto nelem = inpoel.size()/4;
1563 : 2764 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
1564 : 2764 : const auto uncomp = m_u.nprop() / rdof;
1565 : 2764 : const auto pncomp = m_p.nprop() / rdof;
1566 [ + - ]: 5528 : auto ue = m_u;
1567 [ + - ]: 5528 : auto pe = m_p;
1568 : :
1569 : : // If mesh is not refined for field output, cut off ghosts from element
1570 : : // solution. (No need to output ghosts and writer would error.) If mesh is
1571 : : // refined for field output, resize element solution fields to refined mesh.
1572 [ + - ]: 2764 : ue.resize( nelem );
1573 [ + - ]: 2764 : pe.resize( nelem );
1574 : :
1575 : 2764 : auto npoin = coord[0].size();
1576 [ + - ]: 5528 : tk::Fields un( npoin, m_u.nprop()/rdof );
1577 [ + - ]: 5528 : tk::Fields pn( npoin, m_p.nprop()/rdof );
1578 [ + - ]: 2764 : un.fill(0.0);
1579 [ + - ]: 2764 : pn.fill(0.0);
1580 : :
1581 : 2764 : const auto& x = coord[0];
1582 : 2764 : const auto& y = coord[1];
1583 : 2764 : const auto& z = coord[2];
1584 : :
1585 : : // If mesh is not refined for output, evaluate solution in nodes
1586 [ + + ]: 2764 : if (addedTets.empty()) {
1587 : :
1588 [ + + ]: 1330862 : for (std::size_t e=0; e<nelem; ++e) {
1589 : 1328131 : auto e4 = e*4;
1590 : : // Extract element node coordinates
1591 : : std::array< std::array< real, 3>, 4 > ce{{
1592 : 3984393 : {{ x[inpoel[e4 ]], y[inpoel[e4 ]], z[inpoel[e4 ]] }},
1593 : 3984393 : {{ x[inpoel[e4+1]], y[inpoel[e4+1]], z[inpoel[e4+1]] }},
1594 : 3984393 : {{ x[inpoel[e4+2]], y[inpoel[e4+2]], z[inpoel[e4+2]] }},
1595 : 11953179 : {{ x[inpoel[e4+3]], y[inpoel[e4+3]], z[inpoel[e4+3]] }} }};
1596 : : // Compute inverse Jacobian
1597 : 1328131 : auto J = tk::inverseJacobian( ce[0], ce[1], ce[2], ce[3] );
1598 : : // Evaluate solution in child nodes
1599 [ + + ]: 6640655 : for (std::size_t j=0; j<4; ++j) {
1600 : : std::array< real, 3 >
1601 : 5312524 : h{{ce[j][0]-ce[0][0], ce[j][1]-ce[0][1], ce[j][2]-ce[0][2] }};
1602 : 5312524 : auto Bn = tk::eval_basis( m_ndof[e],
1603 [ + - ]: 15937572 : dot(J[0],h), dot(J[1],h), dot(J[2],h) );
1604 [ + - ]: 10625048 : auto u = eval_state( uncomp, 0, rdof, m_ndof[e], e, m_u, Bn, {0, uncomp-1} );
1605 [ + - ]: 10625048 : auto p = eval_state( pncomp, 0, rdof, m_ndof[e], e, m_p, Bn, {0, pncomp-1} );
1606 : : // Assign child node solution
1607 [ + + ][ + - ]: 17638080 : for (std::size_t i=0; i<uncomp; ++i) un(inpoel[e4+j],i,0) += u[i];
1608 [ + + ][ + - ]: 6686076 : for (std::size_t i=0; i<pncomp; ++i) pn(inpoel[e4+j],i,0) += p[i];
1609 : : }
1610 : : }
1611 : :
1612 : : // If mesh is refed for output, evaluate solution in elements and nodes of
1613 : : // refined mesh
1614 : : } else {
1615 : :
1616 [ + - ]: 33 : const auto& pinpoel = Disc()->Inpoel(); // unrefined (parent) mesh
1617 : :
1618 [ + + ]: 157545 : for ([[maybe_unused]] const auto& [child,parent] : addedTets) {
1619 [ - + ][ - - ]: 157512 : Assert( child < nelem, "Indexing out of new solution vector" );
[ - - ][ - - ]
1620 [ - + ][ - - ]: 157512 : Assert( parent < pinpoel.size()/4,
[ - - ][ - - ]
1621 : : "Indexing out of old solution vector" );
1622 : : }
1623 : :
1624 [ + + ]: 157545 : for (const auto& [child,parent] : addedTets) {
1625 : : // Extract parent element's node coordinates
1626 : 157512 : auto p4 = 4*parent;
1627 : : std::array< std::array< real, 3>, 4 > cp{{
1628 : 472536 : {{ x[pinpoel[p4 ]], y[pinpoel[p4 ]], z[pinpoel[p4 ]] }},
1629 : 472536 : {{ x[pinpoel[p4+1]], y[pinpoel[p4+1]], z[pinpoel[p4+1]] }},
1630 : 472536 : {{ x[pinpoel[p4+2]], y[pinpoel[p4+2]], z[pinpoel[p4+2]] }},
1631 : 1417608 : {{ x[pinpoel[p4+3]], y[pinpoel[p4+3]], z[pinpoel[p4+3]] }} }};
1632 : : // Evaluate inverse Jacobian of the parent
1633 : 157512 : auto Jp = tk::inverseJacobian( cp[0], cp[1], cp[2], cp[3] );
1634 : : // Compute child cell centroid
1635 : 157512 : auto c4 = 4*child;
1636 : 157512 : auto cx = (x[inpoel[c4 ]] + x[inpoel[c4+1]] +
1637 : 157512 : x[inpoel[c4+2]] + x[inpoel[c4+3]]) / 4.0;
1638 : 157512 : auto cy = (y[inpoel[c4 ]] + y[inpoel[c4+1]] +
1639 : 157512 : y[inpoel[c4+2]] + y[inpoel[c4+3]]) / 4.0;
1640 : 157512 : auto cz = (z[inpoel[c4 ]] + z[inpoel[c4+1]] +
1641 : 157512 : z[inpoel[c4+2]] + z[inpoel[c4+3]]) / 4.0;
1642 : : // Compute solution in child centroid
1643 : 157512 : std::array< real, 3 > h{{cx-cp[0][0], cy-cp[0][1], cz-cp[0][2] }};
1644 : 157512 : auto B = tk::eval_basis( m_ndof[parent],
1645 [ + - ]: 472536 : dot(Jp[0],h), dot(Jp[1],h), dot(Jp[2],h) );
1646 [ + - ]: 315024 : auto u = eval_state( uncomp, 0, rdof, m_ndof[parent], parent, m_u, B, {0, uncomp-1} );
1647 [ + - ]: 315024 : auto p = eval_state( pncomp, 0, rdof, m_ndof[parent], parent, m_p, B, {0, pncomp-1} );
1648 : : // Assign cell center solution from parent to child
1649 [ + + ][ + - ]: 945072 : for (std::size_t i=0; i<uncomp; ++i) ue(child,i*rdof,0) = u[i];
1650 [ - + ][ - - ]: 157512 : for (std::size_t i=0; i<pncomp; ++i) pe(child,i*rdof,0) = p[i];
1651 : : // Extract child element's node coordinates
1652 : : std::array< std::array< real, 3>, 4 > cc{{
1653 : 472536 : {{ x[inpoel[c4 ]], y[inpoel[c4 ]], z[inpoel[c4 ]] }},
1654 : 472536 : {{ x[inpoel[c4+1]], y[inpoel[c4+1]], z[inpoel[c4+1]] }},
1655 : 472536 : {{ x[inpoel[c4+2]], y[inpoel[c4+2]], z[inpoel[c4+2]] }},
1656 : 1417608 : {{ x[inpoel[c4+3]], y[inpoel[c4+3]], z[inpoel[c4+3]] }} }};
1657 : : // Evaluate solution in child nodes
1658 [ + + ]: 787560 : for (std::size_t j=0; j<4; ++j) {
1659 : : std::array< real, 3 >
1660 : 630048 : hn{{cc[j][0]-cp[0][0], cc[j][1]-cp[0][1], cc[j][2]-cp[0][2] }};
1661 : 630048 : auto Bn = tk::eval_basis( m_ndof[parent],
1662 [ + - ]: 1890144 : dot(Jp[0],hn), dot(Jp[1],hn), dot(Jp[2],hn) );
1663 [ + - ]: 1260096 : auto cnu = eval_state(uncomp, 0, rdof, m_ndof[parent], parent, m_u, Bn, {0, uncomp-1});
1664 [ + - ]: 1260096 : auto cnp = eval_state(pncomp, 0, rdof, m_ndof[parent], parent, m_p, Bn, {0, pncomp-1});
1665 : : // Assign child node solution
1666 [ + + ][ + - ]: 3780288 : for (std::size_t i=0; i<uncomp; ++i) un(inpoel[c4+j],i,0) += cnu[i];
1667 [ - + ][ - - ]: 630048 : for (std::size_t i=0; i<pncomp; ++i) pn(inpoel[c4+j],i,0) += cnp[i];
1668 : : }
1669 : : }
1670 : : }
1671 : :
1672 [ + - ]: 5528 : return { ue, pe, un, pn };
1673 : : }
1674 : :
1675 : : void
1676 : 1026 : DG::lhs()
1677 : : // *****************************************************************************
1678 : : // Compute left-hand side of discrete transport equations
1679 : : // *****************************************************************************
1680 : : {
1681 [ + + ][ + - ]: 2052 : for (const auto& eq : g_dgpde) eq.lhs( m_geoElem, m_lhs );
1682 : :
1683 [ + + ]: 1026 : if (!m_initial) stage();
1684 : 1026 : }
1685 : :
1686 : : void
1687 : 74580 : DG::reco()
1688 : : // *****************************************************************************
1689 : : // Compute reconstructions
1690 : : // *****************************************************************************
1691 : : {
1692 : 74580 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
1693 : 74580 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
1694 : :
1695 : : // Combine own and communicated contributions of unreconstructed solution and
1696 : : // degrees of freedom in cells (if p-adaptive)
1697 [ + + ]: 18892230 : for (const auto& b : m_bid) {
1698 [ - + ][ - - ]: 18817650 : Assert( m_uc[0][b.second].size() == m_u.nprop(), "ncomp size mismatch" );
[ - - ][ - - ]
1699 [ - + ][ - - ]: 18817650 : Assert( m_pc[0][b.second].size() == m_p.nprop(), "ncomp size mismatch" );
[ - - ][ - - ]
1700 [ + + ]: 168233925 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
1701 [ + - ]: 149416275 : m_u(b.first,c,0) = m_uc[0][b.second][c];
1702 : : }
1703 [ + + ]: 38425875 : for (std::size_t c=0; c<m_p.nprop(); ++c) {
1704 [ + - ]: 19608225 : m_p(b.first,c,0) = m_pc[0][b.second][c];
1705 : : }
1706 [ + + ][ + + ]: 18817650 : if (pref && m_stage == 0) {
1707 : 395110 : m_ndof[ b.first ] = m_ndofc[0][ b.second ];
1708 : : }
1709 : : }
1710 : :
1711 [ + + ][ + + ]: 74580 : if (pref && m_stage==0) propagate_ndof();
1712 : :
1713 [ + + ]: 74580 : if (rdof > 1) {
1714 : 44760 : auto d = Disc();
1715 : :
1716 : : // Reconstruct second-order solution and primitive quantities
1717 [ + + ]: 89520 : for (const auto& eq : g_dgpde)
1718 [ + - ]: 44760 : eq.reconstruct( d->T(), m_geoFace, m_geoElem, m_fd, m_esup, m_inpoel,
1719 : 44760 : m_coord, m_u, m_p );
1720 : : }
1721 : :
1722 : : // Send reconstructed solution to neighboring chares
1723 [ + + ]: 74580 : if (m_sendGhost.empty())
1724 : 4395 : comreco_complete();
1725 : : else
1726 [ + + ]: 569715 : for(const auto& [cid, ghostdata] : m_sendGhost) {
1727 [ + - ]: 999060 : std::vector< std::size_t > tetid( ghostdata.size() );
1728 [ + - ]: 999060 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
1729 [ + - ]: 999060 : prim( ghostdata.size() );
1730 : 499530 : std::vector< std::size_t > ndof;
1731 : 499530 : std::size_t j = 0;
1732 [ + + ]: 19317180 : for(const auto& i : ghostdata) {
1733 [ - + ][ - - ]: 18817650 : Assert( i < m_fd.Esuel().size()/4, "Sending reconstructed ghost "
[ - - ][ - - ]
1734 : : "data" );
1735 : 18817650 : tetid[j] = i;
1736 [ + - ]: 18817650 : u[j] = m_u[i];
1737 [ + - ]: 18817650 : prim[j] = m_p[i];
1738 [ + + ][ + + ]: 18817650 : if (pref && m_stage == 0) ndof.push_back( m_ndof[i] );
[ + - ]
1739 : 18817650 : ++j;
1740 : : }
1741 [ + - ][ + - ]: 499530 : thisProxy[ cid ].comreco( thisIndex, tetid, u, prim, ndof );
1742 : : }
1743 : :
1744 : 74580 : ownreco_complete();
1745 : 74580 : }
1746 : :
1747 : : void
1748 : 499530 : DG::comreco( int fromch,
1749 : : const std::vector< std::size_t >& tetid,
1750 : : const std::vector< std::vector< tk::real > >& u,
1751 : : const std::vector< std::vector< tk::real > >& prim,
1752 : : const std::vector< std::size_t >& ndof )
1753 : : // *****************************************************************************
1754 : : // Receive chare-boundary reconstructed ghost data from neighboring chares
1755 : : //! \param[in] fromch Sender chare id
1756 : : //! \param[in] tetid Ghost tet ids we receive solution data for
1757 : : //! \param[in] u Reconstructed high-order solution
1758 : : //! \param[in] prim Limited high-order primitive quantities
1759 : : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
1760 : : //! \details This function receives contributions to the reconstructed solution
1761 : : //! from fellow chares.
1762 : : // *****************************************************************************
1763 : : {
1764 [ - + ][ - - ]: 499530 : Assert( u.size() == tetid.size(), "Size mismatch in DG::comreco()" );
[ - - ][ - - ]
1765 [ - + ][ - - ]: 499530 : Assert( prim.size() == tetid.size(), "Size mismatch in DG::comreco()" );
[ - - ][ - - ]
1766 : :
1767 : 499530 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
1768 : :
1769 [ + + ][ + + ]: 499530 : if (pref && m_stage == 0)
1770 [ - + ][ - - ]: 13280 : Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comreco()" );
[ - - ][ - - ]
1771 : :
1772 : : // Find local-to-ghost tet id map for sender chare
1773 : 499530 : const auto& n = tk::cref_find( m_ghost, fromch );
1774 : :
1775 [ + + ]: 19317180 : for (std::size_t i=0; i<tetid.size(); ++i) {
1776 [ + - ]: 18817650 : auto j = tk::cref_find( n, tetid[i] );
1777 [ - + ][ - - ]: 18817650 : Assert( j >= m_fd.Esuel().size()/4, "Receiving solution non-ghost data" );
[ - - ][ - - ]
1778 [ + - ]: 18817650 : auto b = tk::cref_find( m_bid, j );
1779 [ - + ][ - - ]: 18817650 : Assert( b < m_uc[1].size(), "Indexing out of bounds" );
[ - - ][ - - ]
1780 [ - + ][ - - ]: 18817650 : Assert( b < m_pc[1].size(), "Indexing out of bounds" );
[ - - ][ - - ]
1781 [ + - ]: 18817650 : m_uc[1][b] = u[i];
1782 [ + - ]: 18817650 : m_pc[1][b] = prim[i];
1783 [ + + ][ + + ]: 18817650 : if (pref && m_stage == 0) {
1784 [ - + ][ - - ]: 395110 : Assert( b < m_ndofc[1].size(), "Indexing out of bounds" );
[ - - ][ - - ]
1785 : 395110 : m_ndofc[1][b] = ndof[i];
1786 : : }
1787 : : }
1788 : :
1789 : : // if we have received all solution ghost contributions from neighboring
1790 : : // chares (chares we communicate along chare-boundary faces with), and
1791 : : // contributed our solution to these neighbors, proceed to limiting
1792 [ + + ]: 499530 : if (++m_nreco == m_sendGhost.size()) {
1793 : 70185 : m_nreco = 0;
1794 : 70185 : comreco_complete();
1795 : : }
1796 : 499530 : }
1797 : :
1798 : : void
1799 : 74580 : DG::nodalExtrema()
1800 : : // *****************************************************************************
1801 : : // Compute nodal extrema at chare-boundary nodes. Extrema at internal nodes
1802 : : // are calculated in limiter function.
1803 : : // *****************************************************************************
1804 : : {
1805 [ + - ]: 74580 : auto d = Disc();
1806 [ + - ]: 149160 : auto gid = d->Gid();
1807 [ + - ]: 149160 : auto bid = d->Bid();
1808 : 74580 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
1809 : 74580 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
1810 : 74580 : const auto ncomp = m_u.nprop() / rdof;
1811 : 74580 : const auto nprim = m_p.nprop() / rdof;
1812 : :
1813 : : // Combine own and communicated contributions of unlimited solution, and
1814 : : // if a p-adaptive algorithm is used, degrees of freedom in cells
1815 [ + + ]: 18892230 : for (const auto& [boundary, localtet] : m_bid) {
1816 [ - + ][ - - ]: 18817650 : Assert( m_uc[1][localtet].size() == m_u.nprop(), "ncomp size mismatch" );
[ - - ][ - - ]
1817 [ - + ][ - - ]: 18817650 : Assert( m_pc[1][localtet].size() == m_p.nprop(), "ncomp size mismatch" );
[ - - ][ - - ]
1818 [ + + ]: 168233925 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
1819 [ + - ]: 149416275 : m_u(boundary,c,0) = m_uc[1][localtet][c];
1820 : : }
1821 [ + + ]: 38425875 : for (std::size_t c=0; c<m_p.nprop(); ++c) {
1822 [ + - ]: 19608225 : m_p(boundary,c,0) = m_pc[1][localtet][c];
1823 : : }
1824 [ + + ][ + + ]: 18817650 : if (pref && m_stage == 0) {
1825 : 395110 : m_ndof[ boundary ] = m_ndofc[1][ localtet ];
1826 : : }
1827 : : }
1828 : :
1829 : : // Initialize nodal extrema vector
1830 : 74580 : auto large = std::numeric_limits< tk::real >::max();
1831 [ + + ]: 3905130 : for(std::size_t i = 0; i<bid.size(); i++)
1832 : : {
1833 [ + + ]: 11312280 : for (std::size_t c=0; c<ncomp; ++c)
1834 : : {
1835 [ + + ]: 29926920 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
1836 : : {
1837 : 22445190 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1838 : 22445190 : auto min_mark = max_mark + 1;
1839 : 22445190 : m_uNodalExtrm[i][max_mark] = -large;
1840 : 22445190 : m_uNodalExtrm[i][min_mark] = large;
1841 : : }
1842 : : }
1843 [ + + ]: 4960200 : for (std::size_t c=0; c<nprim; ++c)
1844 : : {
1845 [ + + ]: 4518600 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
1846 : : {
1847 : 3388950 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1848 : 3388950 : auto min_mark = max_mark + 1;
1849 : 3388950 : m_pNodalExtrm[i][max_mark] = -large;
1850 : 3388950 : m_pNodalExtrm[i][min_mark] = large;
1851 : : }
1852 : : }
1853 : : }
1854 : :
1855 : : // Evaluate the max/min value for the chare-boundary nodes
1856 [ + + ]: 74580 : if(rdof > 4) {
1857 [ + - ]: 7785 : evalNodalExtrm(ncomp, nprim, m_ndof_NodalExtrm, d->bndel(), m_inpoel,
1858 [ + - ]: 7785 : m_coord, gid, bid, m_u, m_p, m_uNodalExtrm, m_pNodalExtrm);
1859 : : }
1860 : :
1861 : : // Communicate extrema at nodes to other chares on chare-boundary
1862 [ + + ]: 74580 : if (d->NodeCommMap().empty()) // in serial we are done
1863 [ + - ]: 4395 : comnodalExtrema_complete();
1864 : : else // send nodal extrema to chare-boundary nodes to fellow chares
1865 : : {
1866 [ + + ]: 569715 : for (const auto& [c,n] : d->NodeCommMap()) {
1867 [ + - ][ + - ]: 1498590 : std::vector< std::vector< tk::real > > g1( n.size() ), g2( n.size() );
1868 : 499530 : std::size_t j = 0;
1869 [ + + ]: 6215640 : for (auto i : n)
1870 : : {
1871 [ + - ]: 5716110 : auto p = tk::cref_find(d->Bid(),i);
1872 [ + - ]: 5716110 : g1[ j ] = m_uNodalExtrm[ p ];
1873 [ + - ]: 5716110 : g2[ j++ ] = m_pNodalExtrm[ p ];
1874 : : }
1875 [ + - ][ + - ]: 499530 : thisProxy[c].comnodalExtrema( std::vector<std::size_t>(begin(n),end(n)),
[ + - ]
1876 : : g1, g2 );
1877 : : }
1878 : : }
1879 [ + - ]: 74580 : ownnodalExtrema_complete();
1880 : 74580 : }
1881 : :
1882 : : void
1883 : 499530 : DG::comnodalExtrema( const std::vector< std::size_t >& gid,
1884 : : const std::vector< std::vector< tk::real > >& G1,
1885 : : const std::vector< std::vector< tk::real > >& G2 )
1886 : : // *****************************************************************************
1887 : : // Receive contributions to nodal extrema on chare-boundaries
1888 : : //! \param[in] gid Global mesh node IDs at which we receive grad contributions
1889 : : //! \param[in] G1 Partial contributions of extrema for conservative variables to
1890 : : //! chare-boundary nodes
1891 : : //! \param[in] G2 Partial contributions of extrema for primitive variables to
1892 : : //! chare-boundary nodes
1893 : : //! \details This function receives contributions to m_uNodalExtrm/m_pNodalExtrm
1894 : : //! , which stores nodal extrems at mesh chare-boundary nodes. While
1895 : : //! m_uNodalExtrm/m_pNodalExtrm stores own contributions, m_uNodalExtrmc
1896 : : //! /m_pNodalExtrmc collects the neighbor chare contributions during
1897 : : //! communication.
1898 : : // *****************************************************************************
1899 : : {
1900 [ - + ][ - - ]: 499530 : Assert( G1.size() == gid.size(), "Size mismatch" );
[ - - ][ - - ]
1901 [ - + ][ - - ]: 499530 : Assert( G2.size() == gid.size(), "Size mismatch" );
[ - - ][ - - ]
1902 : :
1903 : 499530 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
1904 : 499530 : const auto ncomp = m_u.nprop() / rdof;
1905 : 499530 : const auto nprim = m_p.nprop() / rdof;
1906 : :
1907 [ + + ]: 6215640 : for (std::size_t i=0; i<gid.size(); ++i)
1908 : : {
1909 : 5716110 : auto& u = m_uNodalExtrmc[gid[i]];
1910 : 5716110 : auto& p = m_pNodalExtrmc[gid[i]];
1911 [ + + ]: 19335360 : for (std::size_t c=0; c<ncomp; ++c)
1912 : : {
1913 [ + + ]: 54477000 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
1914 : : {
1915 : 40857750 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1916 : 40857750 : auto min_mark = max_mark + 1;
1917 : 40857750 : u[max_mark] = std::max( G1[i][max_mark], u[max_mark] );
1918 : 40857750 : u[min_mark] = std::min( G1[i][min_mark], u[min_mark] );
1919 : : }
1920 : : }
1921 [ + + ]: 8418660 : for (std::size_t c=0; c<nprim; ++c)
1922 : : {
1923 [ + + ]: 10810200 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
1924 : : {
1925 : 8107650 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
1926 : 8107650 : auto min_mark = max_mark + 1;
1927 : 8107650 : p[max_mark] = std::max( G2[i][max_mark], p[max_mark] );
1928 : 8107650 : p[min_mark] = std::min( G2[i][min_mark], p[min_mark] );
1929 : : }
1930 : : }
1931 : : }
1932 : :
1933 [ + + ]: 499530 : if (++m_nnodalExtrema == Disc()->NodeCommMap().size())
1934 : : {
1935 : 70185 : m_nnodalExtrema = 0;
1936 : 70185 : comnodalExtrema_complete();
1937 : : }
1938 : 499530 : }
1939 : :
1940 : 75606 : void DG::resizeNodalExtremac()
1941 : : // *****************************************************************************
1942 : : // Resize the buffer vector of nodal extrema
1943 : : // *****************************************************************************
1944 : : {
1945 : 75606 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
1946 : 75606 : const auto ncomp = m_u.nprop() / rdof;
1947 : 75606 : const auto nprim = m_p.nprop() / rdof;
1948 : :
1949 : 75606 : auto large = std::numeric_limits< tk::real >::max();
1950 [ + - ][ + + ]: 583364 : for (const auto& [c,n] : Disc()->NodeCommMap())
1951 : : {
1952 [ + + ]: 6349530 : for (auto i : n) {
1953 [ + - ]: 5841772 : auto& u = m_uNodalExtrmc[i];
1954 [ + - ]: 5841772 : auto& p = m_pNodalExtrmc[i];
1955 [ + - ]: 5841772 : u.resize( 2*m_ndof_NodalExtrm*ncomp, large );
1956 [ + - ]: 5841772 : p.resize( 2*m_ndof_NodalExtrm*nprim, large );
1957 : :
1958 : : // Initialize the minimum nodal extrema
1959 [ + + ]: 23367088 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
1960 : : {
1961 [ + + ]: 59087718 : for(std::size_t k = 0; k < ncomp; k++)
1962 : 41562402 : u[2*k*m_ndof_NodalExtrm+2*idof] = -large;
1963 [ + + ]: 25736712 : for(std::size_t k = 0; k < nprim; k++)
1964 : 8211396 : p[2*k*m_ndof_NodalExtrm+2*idof] = -large;
1965 : : }
1966 : : }
1967 : : }
1968 : 75606 : }
1969 : :
1970 : 7785 : void DG::evalNodalExtrm( const std::size_t ncomp,
1971 : : const std::size_t nprim,
1972 : : const std::size_t ndof_NodalExtrm,
1973 : : const std::vector< std::size_t >& bndel,
1974 : : const std::vector< std::size_t >& inpoel,
1975 : : const tk::UnsMesh::Coords& coord,
1976 : : const std::vector< std::size_t >& gid,
1977 : : const std::unordered_map< std::size_t, std::size_t >&
1978 : : bid,
1979 : : const tk::Fields& U,
1980 : : const tk::Fields& P,
1981 : : std::vector< std::vector<tk::real> >& uNodalExtrm,
1982 : : std::vector< std::vector<tk::real> >& pNodalExtrm )
1983 : : // *****************************************************************************
1984 : : // Compute the nodal extrema for chare-boundary nodes
1985 : : //! \param[in] ncomp Number of conservative variables
1986 : : //! \param[in] nprim Number of primitive variables
1987 : : //! \param[in] ndof_NodalExtrm Degree of freedom for nodal extrema
1988 : : //! \param[in] bndel List of elements contributing to chare-boundary nodes
1989 : : //! \param[in] inpoel Element-node connectivity for element e
1990 : : //! \param[in] coord Array of nodal coordinates
1991 : : //! \param[in] gid Local->global node id map
1992 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
1993 : : //! global node ids (key)
1994 : : //! \param[in] U Vector of conservative variables
1995 : : //! \param[in] P Vector of primitive variables
1996 : : //! \param[in,out] uNodalExtrm Chare-boundary nodal extrema for conservative
1997 : : //! variables
1998 : : //! \param[in,out] pNodalExtrm Chare-boundary nodal extrema for primitive
1999 : : //! variables
2000 : : // *****************************************************************************
2001 : : {
2002 : 7785 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
2003 : :
2004 [ + + ]: 741345 : for (auto e : bndel)
2005 : : {
2006 : : // access node IDs
2007 : : const std::vector<std::size_t> N
2008 [ + - ]: 1467120 : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
2009 : :
2010 : : // Loop over nodes of element e
2011 [ + + ]: 3667800 : for(std::size_t ip=0; ip<4; ++ip)
2012 : : {
2013 [ + - ]: 2934240 : auto i = bid.find( gid[N[ip]] );
2014 [ + + ]: 2934240 : if (i != end(bid)) // If ip is the chare boundary point
2015 : : {
2016 : : // If DG(P2) is applied, find the nodal extrema of the gradients of
2017 : : // conservative/primitive variables in the physical domain
2018 : :
2019 : : // Vector used to store the first order derivatives for both
2020 : : // conservative and primitive variables
2021 [ + - ]: 3893670 : std::vector< std::array< tk::real, 3 > > gradc(ncomp, {0.0, 0.0, 0.0});
2022 [ + - ]: 3893670 : std::vector< std::array< tk::real, 3 > > gradp(ncomp, {0.0, 0.0, 0.0});
2023 : :
2024 : 1946835 : const auto& cx = coord[0];
2025 : 1946835 : const auto& cy = coord[1];
2026 : 1946835 : const auto& cz = coord[2];
2027 : :
2028 : : std::array< std::array< tk::real, 3>, 4 > coordel {{
2029 : 5840505 : {{ cx[ N[0] ], cy[ N[0] ], cz[ N[0] ] }},
2030 : 5840505 : {{ cx[ N[1] ], cy[ N[1] ], cz[ N[1] ] }},
2031 : 5840505 : {{ cx[ N[2] ], cy[ N[2] ], cz[ N[2] ] }},
2032 : 5840505 : {{ cx[ N[3] ], cy[ N[3] ], cz[ N[3] ] }}
2033 : 21415185 : }};
2034 : :
2035 : 1946835 : auto jacInv = tk::inverseJacobian( coordel[0], coordel[1],
2036 : 3893670 : coordel[2], coordel[3] );
2037 : :
2038 : : // Compute the derivatives of basis functions
2039 [ + - ]: 3893670 : auto dBdx = tk::eval_dBdx_p1( rdof, jacInv );
2040 : :
2041 : 3893670 : std::array< std::vector< tk::real >, 3 > center;
2042 [ + - ]: 1946835 : center[0].resize(1, 0.25);
2043 [ + - ]: 1946835 : center[1].resize(1, 0.25);
2044 [ + - ]: 1946835 : center[2].resize(1, 0.25);
2045 [ + - ]: 1946835 : tk::eval_dBdx_p2(0, center, jacInv, dBdx);
2046 : :
2047 : : // Evaluate the first order derivative in physical domain
2048 [ + + ]: 11254110 : for(std::size_t icomp = 0; icomp < ncomp; icomp++)
2049 : : {
2050 : 9307275 : auto mark = icomp * rdof;
2051 [ + + ]: 37229100 : for(std::size_t idir = 0; idir < 3; idir++)
2052 : : {
2053 : 27921825 : gradc[icomp][idir] = 0;
2054 [ + + ]: 279218250 : for(std::size_t idof = 1; idof < rdof; idof++)
2055 [ + - ]: 251296425 : gradc[icomp][idir] += U(e, mark+idof, 0) * dBdx[idir][idof];
2056 : : }
2057 : : }
2058 [ - + ]: 1946835 : for(std::size_t icomp = 0; icomp < nprim; icomp++)
2059 : : {
2060 : 0 : auto mark = icomp * rdof;
2061 [ - - ]: 0 : for(std::size_t idir = 0; idir < 3; idir++)
2062 : : {
2063 : 0 : gradp[icomp][idir] = 0;
2064 [ - - ]: 0 : for(std::size_t idof = 1; idof < rdof; idof++)
2065 [ - - ]: 0 : gradp[icomp][idir] += P(e, mark+idof, 0) * dBdx[idir][idof];
2066 : : }
2067 : : }
2068 : :
2069 : : // Store the extrema for the gradients
2070 [ + + ]: 11254110 : for (std::size_t c=0; c<ncomp; ++c)
2071 : : {
2072 [ + + ]: 37229100 : for (std::size_t idof = 0; idof < ndof_NodalExtrm; idof++)
2073 : : {
2074 : 27921825 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
2075 : 27921825 : auto min_mark = max_mark + 1;
2076 : 27921825 : auto& ex = uNodalExtrm[i->second];
2077 : 27921825 : ex[max_mark] = std::max(ex[max_mark], gradc[c][idof-1]);
2078 : 27921825 : ex[min_mark] = std::min(ex[min_mark], gradc[c][idof-1]);
2079 : : }
2080 : : }
2081 [ - + ]: 1946835 : for (std::size_t c=0; c<nprim; ++c)
2082 : : {
2083 [ - - ]: 0 : for (std::size_t idof = 0; idof < ndof_NodalExtrm; idof++)
2084 : : {
2085 : 0 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
2086 : 0 : auto min_mark = max_mark + 1;
2087 : 0 : auto& ex = pNodalExtrm[i->second];
2088 : 0 : ex[max_mark] = std::max(ex[max_mark], gradp[c][idof-1]);
2089 : 0 : ex[min_mark] = std::min(ex[min_mark], gradp[c][idof-1]);
2090 : : }
2091 : : }
2092 : : }
2093 : : }
2094 : : }
2095 : 7785 : }
2096 : :
2097 : : void
2098 : 74580 : DG::lim()
2099 : : // *****************************************************************************
2100 : : // Compute limiter function
2101 : : // *****************************************************************************
2102 : : {
2103 : 74580 : auto d = Disc();
2104 : 74580 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
2105 : 74580 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
2106 : 74580 : const auto ncomp = m_u.nprop() / rdof;
2107 : 74580 : const auto nprim = m_p.nprop() / rdof;
2108 : :
2109 : : // Combine own and communicated contributions to nodal extrema
2110 [ + + ]: 3905130 : for (const auto& [gid,g] : m_uNodalExtrmc) {
2111 [ + - ]: 3830550 : auto bid = tk::cref_find( d->Bid(), gid );
2112 [ + + ]: 11312280 : for (ncomp_t c=0; c<ncomp; ++c)
2113 : : {
2114 [ + + ]: 29926920 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
2115 : : {
2116 : 22445190 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
2117 : 22445190 : auto min_mark = max_mark + 1;
2118 : 22445190 : m_uNodalExtrm[bid][max_mark] =
2119 : 22445190 : std::max(g[max_mark], m_uNodalExtrm[bid][max_mark]);
2120 : 22445190 : m_uNodalExtrm[bid][min_mark] =
2121 : 22445190 : std::min(g[min_mark], m_uNodalExtrm[bid][min_mark]);
2122 : : }
2123 : : }
2124 : : }
2125 [ + + ]: 3905130 : for (const auto& [gid,g] : m_pNodalExtrmc) {
2126 [ + - ]: 3830550 : auto bid = tk::cref_find( d->Bid(), gid );
2127 [ + + ]: 4960200 : for (ncomp_t c=0; c<nprim; ++c)
2128 : : {
2129 [ + + ]: 4518600 : for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
2130 : : {
2131 : 3388950 : auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
2132 : 3388950 : auto min_mark = max_mark + 1;
2133 : 3388950 : m_pNodalExtrm[bid][max_mark] =
2134 : 3388950 : std::max(g[max_mark], m_pNodalExtrm[bid][max_mark]);
2135 : 3388950 : m_pNodalExtrm[bid][min_mark] =
2136 : 3388950 : std::min(g[min_mark], m_pNodalExtrm[bid][min_mark]);
2137 : : }
2138 : : }
2139 : : }
2140 : :
2141 : : // clear gradients receive buffer
2142 : 74580 : tk::destroy(m_uNodalExtrmc);
2143 : 74580 : tk::destroy(m_pNodalExtrmc);
2144 : :
2145 [ + + ]: 74580 : if (rdof > 1)
2146 [ + + ]: 89520 : for (const auto& eq : g_dgpde)
2147 [ + - ]: 89520 : eq.limit( d->T(), m_geoFace, m_geoElem, m_fd, m_esup, m_inpoel, m_coord,
2148 : 44760 : m_ndof, d->Gid(), d->Bid(), m_uNodalExtrm, m_pNodalExtrm, m_u,
2149 : 44760 : m_p, m_shockmarker );
2150 : :
2151 : : // Send limited solution to neighboring chares
2152 [ + + ]: 74580 : if (m_sendGhost.empty())
2153 : 4395 : comlim_complete();
2154 : : else
2155 [ + + ]: 569715 : for(const auto& [cid, ghostdata] : m_sendGhost) {
2156 [ + - ]: 999060 : std::vector< std::size_t > tetid( ghostdata.size() );
2157 [ + - ]: 999060 : std::vector< std::vector< tk::real > > u( ghostdata.size() ),
2158 [ + - ]: 999060 : prim( ghostdata.size() );
2159 : 499530 : std::vector< std::size_t > ndof;
2160 : 499530 : std::size_t j = 0;
2161 [ + + ]: 19317180 : for(const auto& i : ghostdata) {
2162 [ - + ][ - - ]: 18817650 : Assert( i < m_fd.Esuel().size()/4, "Sending limiter ghost data" );
[ - - ][ - - ]
2163 : 18817650 : tetid[j] = i;
2164 [ + - ]: 18817650 : u[j] = m_u[i];
2165 [ + - ]: 18817650 : prim[j] = m_p[i];
2166 [ + + ][ + + ]: 18817650 : if (pref && m_stage == 0) ndof.push_back( m_ndof[i] );
[ + - ]
2167 : 18817650 : ++j;
2168 : : }
2169 [ + - ][ + - ]: 499530 : thisProxy[ cid ].comlim( thisIndex, tetid, u, prim, ndof );
2170 : : }
2171 : :
2172 : 74580 : ownlim_complete();
2173 : 74580 : }
2174 : :
2175 : : void
2176 : 1770 : DG::propagate_ndof()
2177 : : // *****************************************************************************
2178 : : // p-refine all elements that are adjacent to p-refined elements
2179 : : //! \details This function p-refines all the neighbors of an element that has
2180 : : //! been p-refined as a result of an error indicator.
2181 : : // *****************************************************************************
2182 : : {
2183 : 1770 : const auto& esuf = m_fd.Esuf();
2184 : :
2185 : : // Copy number of degrees of freedom for each cell
2186 [ + - ]: 3540 : auto ndof = m_ndof;
2187 : :
2188 : : // p-refine all neighboring elements of elements that have been p-refined as a
2189 : : // result of error indicators
2190 [ + - ][ + + ]: 753210 : for( auto f=m_fd.Nbfac(); f<esuf.size()/2; ++f )
2191 : : {
2192 : 751440 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
2193 : 751440 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
2194 : :
2195 [ + + ]: 751440 : if (m_ndof[el] > m_ndof[er])
2196 : 5111 : ndof[er] = m_ndof[el];
2197 : :
2198 [ + + ]: 751440 : if (m_ndof[el] < m_ndof[er])
2199 : 4832 : ndof[el] = m_ndof[er];
2200 : : }
2201 : :
2202 : : // Update number of degrees of freedom for each cell
2203 [ + - ]: 1770 : m_ndof = ndof;
2204 : 1770 : }
2205 : :
2206 : : void
2207 : 499530 : DG::comlim( int fromch,
2208 : : const std::vector< std::size_t >& tetid,
2209 : : const std::vector< std::vector< tk::real > >& u,
2210 : : const std::vector< std::vector< tk::real > >& prim,
2211 : : const std::vector< std::size_t >& ndof )
2212 : : // *****************************************************************************
2213 : : // Receive chare-boundary limiter ghost data from neighboring chares
2214 : : //! \param[in] fromch Sender chare id
2215 : : //! \param[in] tetid Ghost tet ids we receive solution data for
2216 : : //! \param[in] u Limited high-order solution
2217 : : //! \param[in] prim Limited high-order primitive quantities
2218 : : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
2219 : : //! \details This function receives contributions to the limited solution from
2220 : : //! fellow chares.
2221 : : // *****************************************************************************
2222 : : {
2223 [ - + ][ - - ]: 499530 : Assert( u.size() == tetid.size(), "Size mismatch in DG::comlim()" );
[ - - ][ - - ]
2224 [ - + ][ - - ]: 499530 : Assert( prim.size() == tetid.size(), "Size mismatch in DG::comlim()" );
[ - - ][ - - ]
2225 : :
2226 : 499530 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
2227 : :
2228 [ + + ][ + + ]: 499530 : if (pref && m_stage == 0)
2229 [ - + ][ - - ]: 13280 : Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comlim()" );
[ - - ][ - - ]
2230 : :
2231 : : // Find local-to-ghost tet id map for sender chare
2232 : 499530 : const auto& n = tk::cref_find( m_ghost, fromch );
2233 : :
2234 [ + + ]: 19317180 : for (std::size_t i=0; i<tetid.size(); ++i) {
2235 [ + - ]: 18817650 : auto j = tk::cref_find( n, tetid[i] );
2236 [ - + ][ - - ]: 18817650 : Assert( j >= m_fd.Esuel().size()/4, "Receiving solution non-ghost data" );
[ - - ][ - - ]
2237 [ + - ]: 18817650 : auto b = tk::cref_find( m_bid, j );
2238 [ - + ][ - - ]: 18817650 : Assert( b < m_uc[2].size(), "Indexing out of bounds" );
[ - - ][ - - ]
2239 [ - + ][ - - ]: 18817650 : Assert( b < m_pc[2].size(), "Indexing out of bounds" );
[ - - ][ - - ]
2240 [ + - ]: 18817650 : m_uc[2][b] = u[i];
2241 [ + - ]: 18817650 : m_pc[2][b] = prim[i];
2242 [ + + ][ + + ]: 18817650 : if (pref && m_stage == 0) {
2243 [ - + ][ - - ]: 395110 : Assert( b < m_ndofc[2].size(), "Indexing out of bounds" );
[ - - ][ - - ]
2244 : 395110 : m_ndofc[2][b] = ndof[i];
2245 : : }
2246 : : }
2247 : :
2248 : : // if we have received all solution ghost contributions from neighboring
2249 : : // chares (chares we communicate along chare-boundary faces with), and
2250 : : // contributed our solution to these neighbors, proceed to limiting
2251 [ + + ]: 499530 : if (++m_nlim == m_sendGhost.size()) {
2252 : 70185 : m_nlim = 0;
2253 : 70185 : comlim_complete();
2254 : : }
2255 : 499530 : }
2256 : :
2257 : : void
2258 : 74580 : DG::dt()
2259 : : // *****************************************************************************
2260 : : // Compute time step size
2261 : : // *****************************************************************************
2262 : : {
2263 : 74580 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
2264 [ + - ]: 74580 : auto d = Disc();
2265 : :
2266 : :
2267 : : // Combine own and communicated contributions of limited solution and degrees
2268 : : // of freedom in cells (if p-adaptive)
2269 [ + + ]: 18892230 : for (const auto& b : m_bid) {
2270 [ - + ][ - - ]: 18817650 : Assert( m_uc[2][b.second].size() == m_u.nprop(), "ncomp size mismatch" );
[ - - ][ - - ]
2271 [ - + ][ - - ]: 18817650 : Assert( m_pc[2][b.second].size() == m_p.nprop(), "ncomp size mismatch" );
[ - - ][ - - ]
2272 [ + + ]: 168233925 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
2273 [ + - ]: 149416275 : m_u(b.first,c,0) = m_uc[2][b.second][c];
2274 : : }
2275 [ + + ]: 38425875 : for (std::size_t c=0; c<m_p.nprop(); ++c) {
2276 [ + - ]: 19608225 : m_p(b.first,c,0) = m_pc[2][b.second][c];
2277 : : }
2278 [ + + ][ + + ]: 18817650 : if (pref && m_stage == 0) {
2279 : 395110 : m_ndof[ b.first ] = m_ndofc[2][ b.second ];
2280 : : }
2281 : : }
2282 : :
2283 : 74580 : auto mindt = std::numeric_limits< tk::real >::max();
2284 : :
2285 [ + + ]: 74580 : if (m_stage == 0)
2286 : : {
2287 : 24860 : auto const_dt = g_inputdeck.get< tag::discr, tag::dt >();
2288 : 24860 : auto def_const_dt = g_inputdeck_defaults.get< tag::discr, tag::dt >();
2289 : 24860 : auto eps = std::numeric_limits< tk::real >::epsilon();
2290 : :
2291 : : // use constant dt if configured
2292 [ + + ]: 24860 : if (std::abs(const_dt - def_const_dt) > eps) {
2293 : :
2294 : 22415 : mindt = const_dt;
2295 : :
2296 : : } else { // compute dt based on CFL
2297 : :
2298 : : // find the minimum dt across all PDEs integrated
2299 [ + + ]: 4890 : for (const auto& eq : g_dgpde) {
2300 : : auto eqdt =
2301 : 4890 : eq.dt( m_coord, m_inpoel, m_fd, m_geoFace, m_geoElem, m_ndof,
2302 [ + - ]: 2445 : m_u, m_p, m_fd.Esuel().size()/4 );
2303 [ + - ]: 2445 : if (eqdt < mindt) mindt = eqdt;
2304 : : }
2305 : :
2306 : 2445 : mindt *= g_inputdeck.get< tag::discr, tag::cfl >();
2307 : : }
2308 : : }
2309 : : else
2310 : : {
2311 : 49720 : mindt = d->Dt();
2312 : : }
2313 : :
2314 : : // Resize the buffer vector of nodal extrema
2315 [ + - ]: 74580 : resizeNodalExtremac();
2316 : :
2317 : : // Contribute to minimum dt across all chares then advance to next step
2318 [ + - ]: 74580 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
2319 [ + - ][ + - ]: 149160 : CkCallback(CkReductionTarget(DG,solve), thisProxy) );
2320 : 74580 : }
2321 : :
2322 : : void
2323 : 74580 : DG::solve( tk::real newdt )
2324 : : // *****************************************************************************
2325 : : // Compute right-hand side of discrete transport equations
2326 : : //! \param[in] newdt Size of this new time step
2327 : : // *****************************************************************************
2328 : : {
2329 : : // Enable SDAG wait for building the solution vector during the next stage
2330 [ + - ]: 74580 : thisProxy[ thisIndex ].wait4sol();
2331 [ + - ]: 74580 : thisProxy[ thisIndex ].wait4reco();
2332 [ + - ]: 74580 : thisProxy[ thisIndex ].wait4nodalExtrema();
2333 [ + - ]: 74580 : thisProxy[ thisIndex ].wait4lim();
2334 [ + - ]: 74580 : thisProxy[ thisIndex ].wait4nod();
2335 : :
2336 : 74580 : auto d = Disc();
2337 : 74580 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
2338 : 74580 : const auto ndof = g_inputdeck.get< tag::discr, tag::ndof >();
2339 : 74580 : const auto neq = m_u.nprop()/rdof;
2340 : :
2341 : : // Set new time step size
2342 [ + + ]: 74580 : if (m_stage == 0) d->setdt( newdt );
2343 : :
2344 : 74580 : const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
2345 [ + + ][ + + ]: 74580 : if (pref && m_stage == 0)
2346 : : {
2347 : : // When the element are coarsened, high order terms should be zero
2348 [ + + ]: 809600 : for(std::size_t e = 0; e < m_nunk; e++)
2349 : : {
2350 : 807830 : const auto ncomp= m_u.nprop()/rdof;
2351 [ + + ]: 807830 : if(m_ndof[e] == 1)
2352 : : {
2353 [ + + ]: 3824190 : for (std::size_t c=0; c<ncomp; ++c)
2354 : : {
2355 : 3186825 : auto mark = c*rdof;
2356 : 3186825 : m_u(e, mark+1, 0) = 0.0;
2357 : 3186825 : m_u(e, mark+2, 0) = 0.0;
2358 : 3186825 : m_u(e, mark+3, 0) = 0.0;
2359 : : }
2360 [ + + ]: 170465 : } else if(m_ndof[e] == 4)
2361 : : {
2362 [ + + ]: 498132 : for (std::size_t c=0; c<ncomp; ++c)
2363 : : {
2364 : 415110 : auto mark = c*ndof;
2365 : 415110 : m_u(e, mark+4, 0) = 0.0;
2366 : 415110 : m_u(e, mark+5, 0) = 0.0;
2367 : 415110 : m_u(e, mark+6, 0) = 0.0;
2368 : 415110 : m_u(e, mark+7, 0) = 0.0;
2369 : 415110 : m_u(e, mark+8, 0) = 0.0;
2370 : 415110 : m_u(e, mark+9, 0) = 0.0;
2371 : : }
2372 : : }
2373 : : }
2374 : : }
2375 : :
2376 : : // Update Un
2377 [ + + ]: 74580 : if (m_stage == 0) m_un = m_u;
2378 : :
2379 [ + + ]: 149160 : for (const auto& eq : g_dgpde)
2380 [ + - ]: 74580 : eq.rhs( d->T(), m_geoFace, m_geoElem, m_fd, m_inpoel, m_boxelems, m_coord,
2381 : 74580 : m_u, m_p, m_ndof, m_rhs );
2382 : :
2383 : : // Explicit time-stepping using RK3 to discretize time-derivative
2384 [ + + ]: 54466800 : for(std::size_t e=0; e<m_nunk; ++e)
2385 [ + + ]: 170108760 : for(std::size_t c=0; c<neq; ++c)
2386 : : {
2387 [ + + ]: 457366830 : for (std::size_t k=0; k<m_numEqDof[c]; ++k)
2388 : : {
2389 : 341650290 : auto rmark = c*rdof+k;
2390 : 341650290 : auto mark = c*ndof+k;
2391 : 341650290 : m_u(e, rmark, 0) = rkcoef[0][m_stage] * m_un(e, rmark, 0)
2392 : 341650290 : + rkcoef[1][m_stage] * ( m_u(e, rmark, 0)
2393 : 341650290 : + d->Dt() * m_rhs(e, mark, 0)/m_lhs(e, mark, 0) );
2394 [ + + ]: 341650290 : if(fabs(m_u(e, rmark, 0)) < 1e-16)
2395 : 170241950 : m_u(e, rmark, 0) = 0;
2396 : : }
2397 : : // zero out unused/reconstructed dofs of equations using reduced dofs
2398 : : // (see DGMultiMat::numEquationDofs())
2399 [ + + ]: 115716540 : if (m_numEqDof[c] < rdof) {
2400 [ + + ]: 66551100 : for (std::size_t k=m_numEqDof[c]; k<rdof; ++k)
2401 : : {
2402 : 49913325 : auto rmark = c*rdof+k;
2403 : 49913325 : m_u(e, rmark, 0) = 0.0;
2404 : : }
2405 : : }
2406 : : }
2407 : :
2408 : : // Update primitives based on the evolved solution
2409 [ + + ]: 149160 : for (const auto& eq : g_dgpde)
2410 : : {
2411 [ + - ]: 74580 : eq.updateInterfaceCells( m_u, m_fd.Esuel().size()/4, m_ndof );
2412 [ + - ]: 74580 : eq.updatePrimitives( m_u, m_lhs, m_geoElem, m_p, m_fd.Esuel().size()/4 );
2413 [ + - ]: 74580 : eq.cleanTraceMaterial( m_geoElem, m_u, m_p, m_fd.Esuel().size()/4 );
2414 : : }
2415 : :
2416 [ + + ]: 74580 : if (m_stage < 2) {
2417 : :
2418 : : // continue with next time step stage
2419 : 49720 : stage();
2420 : :
2421 : : } else {
2422 : :
2423 : : // Compute diagnostics, e.g., residuals
2424 : 24860 : auto diag_computed = m_diag.compute( *d, m_u.nunk()-m_fd.Esuel().size()/4,
2425 : 24860 : m_geoElem, m_ndof, m_u );
2426 : :
2427 : : // Increase number of iterations and physical time
2428 : 24860 : d->next();
2429 : :
2430 : : // Continue to mesh refinement (if configured)
2431 [ + + ][ + - ]: 24860 : if (!diag_computed) refine( std::vector< tk::real >( m_u.nprop(), 0.0 ) );
[ + - ]
2432 : :
2433 : : }
2434 : 74580 : }
2435 : :
2436 : : void
2437 : 24860 : DG::refine( [[maybe_unused]] const std::vector< tk::real >& l2res )
2438 : : // *****************************************************************************
2439 : : // Optionally refine/derefine mesh
2440 : : //! \param[in] l2res L2-norms of the residual for each scalar component
2441 : : //! computed across the whole problem
2442 : : // *****************************************************************************
2443 : : {
2444 : 24860 : auto d = Disc();
2445 : :
2446 : 24860 : auto dtref = g_inputdeck.get< tag::amr, tag::dtref >();
2447 : 24860 : auto dtfreq = g_inputdeck.get< tag::amr, tag::dtfreq >();
2448 : :
2449 : : // if t>0 refinement enabled and we hit the dtref frequency
2450 [ + + ][ + + ]: 24860 : if (dtref && !(d->It() % dtfreq)) { // refine
[ + + ]
2451 : :
2452 : 108 : d->startvol();
2453 [ + - ]: 108 : d->Ref()->dtref( m_fd.Bface(), {}, tk::remap(m_fd.Triinpoel(),d->Gid()) );
2454 : 108 : d->refined() = 1;
2455 : :
2456 : : } else { // do not refine
2457 : :
2458 : 24752 : d->refined() = 0;
2459 : 24752 : stage();
2460 : :
2461 : : }
2462 : 24860 : }
2463 : :
2464 : : void
2465 : 108 : DG::resizePostAMR(
2466 : : const std::vector< std::size_t >& /*ginpoel*/,
2467 : : const tk::UnsMesh::Chunk& chunk,
2468 : : const tk::UnsMesh::Coords& coord,
2469 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /*addedNodes*/,
2470 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
2471 : : const std::set< std::size_t >& /*removedNodes*/,
2472 : : const tk::NodeCommMap& nodeCommMap,
2473 : : const std::map< int, std::vector< std::size_t > >& bface,
2474 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
2475 : : const std::vector< std::size_t >& triinpoel )
2476 : : // *****************************************************************************
2477 : : // Receive new mesh from Refiner
2478 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
2479 : : //! \param[in] coord New mesh node coordinates
2480 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
2481 : : //! \param[in] nodeCommMap New node communication map
2482 : : //! \param[in] bface Boundary-faces mapped to side set ids
2483 : : //! \param[in] triinpoel Boundary-face connectivity
2484 : : // *****************************************************************************
2485 : : {
2486 [ + - ]: 108 : auto d = Disc();
2487 : :
2488 : : // Set flag that indicates that we are during time stepping
2489 : 108 : m_initial = 0;
2490 : :
2491 : : // Zero field output iteration count between two mesh refinement steps
2492 : 108 : d->Itf() = 0;
2493 : :
2494 : : // Increase number of iterations with mesh refinement
2495 : 108 : ++d->Itr();
2496 : :
2497 : : // Save old number of elements
2498 : 108 : [[maybe_unused]] auto old_nelem = m_inpoel.size()/4;
2499 : :
2500 : : // Resize mesh data structures
2501 [ + - ]: 108 : d->resizePostAMR( chunk, coord, nodeCommMap );
2502 : :
2503 : : // Update state
2504 [ + - ]: 108 : m_inpoel = d->Inpoel();
2505 [ + - ]: 108 : m_coord = d->Coord();
2506 : 108 : auto nelem = m_inpoel.size()/4;
2507 [ + - ]: 108 : m_p.resize( nelem );
2508 [ + - ]: 108 : m_u.resize( nelem );
2509 [ + - ]: 108 : m_un.resize( nelem );
2510 [ + - ]: 108 : m_lhs.resize( nelem );
2511 [ + - ]: 108 : m_rhs.resize( nelem );
2512 [ + - ][ + - ]: 216 : m_uNodalExtrm.resize( Disc()->Bid().size(), std::vector<tk::real>( 2*
2513 [ + - ]: 108 : m_ndof_NodalExtrm*g_inputdeck.get< tag::component >().nprop() ) );
2514 [ + - ][ + - ]: 216 : m_pNodalExtrm.resize( Disc()->Bid().size(), std::vector<tk::real>( 2*
2515 [ + - ]: 108 : m_ndof_NodalExtrm*m_p.nprop()/g_inputdeck.get< tag::discr, tag::rdof >()));
2516 : :
2517 : : // Resize the buffer vector of nodal extrema
2518 [ + - ]: 108 : resizeNodalExtremac();
2519 : :
2520 [ + - ][ + - ]: 108 : m_fd = FaceData( m_inpoel, bface, tk::remap(triinpoel,d->Lid()) );
2521 : :
2522 : : m_geoFace =
2523 [ + - ]: 108 : tk::Fields( tk::genGeoFaceTri( m_fd.Nipfac(), m_fd.Inpofa(), coord ) );
2524 [ + - ]: 108 : m_geoElem = tk::Fields( tk::genGeoElemTet( m_inpoel, coord ) );
2525 : :
2526 : 108 : m_nfac = m_fd.Inpofa().size()/3;
2527 : 108 : m_nunk = nelem;
2528 : 108 : m_npoin = coord[0].size();
2529 : 108 : m_bndFace.clear();
2530 : 108 : m_exptGhost.clear();
2531 : 108 : m_sendGhost.clear();
2532 : 108 : m_ghost.clear();
2533 : 108 : m_esup.clear();
2534 : :
2535 : : // Update solution on new mesh, P0 (cell center value) only for now
2536 [ + - ]: 108 : m_un = m_u;
2537 [ + - ]: 216 : auto pn = m_p;
2538 : 108 : auto unprop = m_u.nprop();
2539 : 108 : auto pnprop = m_p.nprop();
2540 [ + + ]: 64620 : for (const auto& [child,parent] : addedTets) {
2541 [ - + ][ - - ]: 64512 : Assert( child < nelem, "Indexing out of new solution vector" );
[ - - ][ - - ]
2542 [ - + ][ - - ]: 64512 : Assert( parent < old_nelem, "Indexing out of old solution vector" );
[ - - ][ - - ]
2543 [ + + ][ + - ]: 129024 : for (std::size_t i=0; i<unprop; ++i) m_u(child,i,0) = m_un(parent,i,0);
[ + - ]
2544 [ - + ][ - - ]: 64512 : for (std::size_t i=0; i<pnprop; ++i) m_p(child,i,0) = pn(parent,i,0);
[ - - ]
2545 : : }
2546 [ + - ]: 108 : m_un = m_u;
2547 : :
2548 : : // Enable SDAG wait for setting up chare boundary faces
2549 [ + - ][ + - ]: 108 : thisProxy[ thisIndex ].wait4fac();
2550 : :
2551 : : // Resize communication buffers
2552 [ + - ]: 108 : resizeComm();
2553 : 108 : }
2554 : :
2555 : : bool
2556 : 22434 : DG::fieldOutput() const
2557 : : // *****************************************************************************
2558 : : // Decide wether to output field data
2559 : : //! \return True if field data is output in this step
2560 : : // *****************************************************************************
2561 : : {
2562 : 22434 : auto d = Disc();
2563 : :
2564 : : // Output field data
2565 [ + + ][ + - ]: 22434 : return d->fielditer() or d->fieldtime() or d->fieldrange() or d->finished();
[ + - ][ - + ]
2566 : : }
2567 : :
2568 : : bool
2569 : 2764 : DG::refinedOutput() const
2570 : : // *****************************************************************************
2571 : : // Decide if we write field output using a refined mesh
2572 : : //! \return True if field output will use a refined mesh
2573 : : // *****************************************************************************
2574 : : {
2575 [ + + ]: 2797 : return g_inputdeck.get< tag::cmd, tag::io, tag::refined >() &&
2576 [ + - ]: 2797 : g_inputdeck.get< tag::discr, tag::scheme >() != ctr::SchemeType::DG;
2577 : : }
2578 : :
2579 : : void
2580 : 2764 : DG::writeFields( CkCallback c )
2581 : : // *****************************************************************************
2582 : : // Output mesh field data
2583 : : //! \param[in] c Function to continue with after the write
2584 : : // *****************************************************************************
2585 : : {
2586 [ + - ]: 2764 : auto d = Disc();
2587 : :
2588 : : // Output time history
2589 [ + - ][ + - ]: 2764 : if (d->histiter() or d->histtime() or d->histrange()) {
[ + - ][ + - ]
[ + - ][ - + ]
[ - + ]
2590 : 0 : std::vector< std::vector< tk::real > > hist;
2591 [ - - ]: 0 : for (const auto& eq : g_dgpde) {
2592 [ - - ]: 0 : auto h = eq.histOutput( d->Hist(), m_inpoel, m_coord, m_u, m_p );
2593 [ - - ]: 0 : hist.insert( end(hist), begin(h), end(h) );
2594 : : }
2595 [ - - ]: 0 : d->history( std::move(hist) );
2596 : : }
2597 : :
2598 : 2764 : const auto& inpoel = std::get< 0 >( m_outmesh.chunk );
2599 [ + - ]: 5528 : auto esup = tk::genEsup( inpoel, 4 );
2600 : :
2601 : : // Combine own and communicated contributions and finish averaging of node
2602 : : // field output in chare boundary nodes
2603 : 2764 : const auto& lid = std::get< 2 >( m_outmesh.chunk );
2604 [ + + ]: 105140 : for (const auto& [g,f] : m_nodefieldsc) {
2605 [ - + ][ - - ]: 102376 : Assert( m_nodefields.size() == f.first.size(), "Size mismatch" );
[ - - ][ - - ]
2606 [ + - ]: 102376 : auto p = tk::cref_find( lid, g );
2607 [ - + ]: 102376 : for (std::size_t i=0; i<f.first.size(); ++i) {
2608 : 0 : m_nodefields[i][p] += f.first[i];
2609 : 0 : m_nodefields[i][p] /= static_cast< tk::real >(
2610 : 0 : esup.second[p+1] - esup.second[p] + f.second );
2611 : : }
2612 : : }
2613 : 2764 : tk::destroy( m_nodefieldsc );
2614 : :
2615 : : // Lambda to decide if a node (global id) is on a chare boundary of the field
2616 : : // output mesh. p - global node id, return true if node is on the chare
2617 : : // boundary.
2618 : 898128 : auto chbnd = [ this ]( std::size_t p ) {
2619 : : return
2620 : 449064 : std::any_of( m_outmesh.nodeCommMap.cbegin(), m_outmesh.nodeCommMap.cend(),
2621 [ + - ]: 1276352 : [&](const auto& s) { return s.second.find(p) != s.second.cend(); } );
2622 : 2764 : };
2623 : :
2624 : : // Finish computing node field output averages in internal nodes
2625 : 2764 : auto npoin = m_outmesh.coord[0].size();
2626 : 2764 : auto& gid = std::get< 1 >( m_outmesh.chunk );
2627 [ + + ]: 451828 : for (std::size_t p=0; p<npoin; ++p) {
2628 [ + - ][ + + ]: 449064 : if (!chbnd(gid[p])) {
2629 : 346688 : auto n = static_cast< tk::real >( esup.second[p+1] - esup.second[p] );
2630 [ - + ]: 346688 : for (auto& f : m_nodefields) f[p] /= n;
2631 : : }
2632 : : }
2633 : :
2634 : : // Query fields names requested by user
2635 [ + - ]: 5528 : auto elemfieldnames = numericFieldNames( tk::Centering::ELEM );
2636 [ + - ]: 2764 : auto nodefieldnames = numericFieldNames( tk::Centering::NODE );
2637 : :
2638 : : // Collect field output names for analytical solutions
2639 [ + + ]: 5528 : for (const auto& eq : g_dgpde) {
2640 [ + - ]: 2764 : analyticFieldNames( eq, tk::Centering::ELEM, elemfieldnames );
2641 [ + - ]: 2764 : analyticFieldNames( eq, tk::Centering::NODE, nodefieldnames );
2642 : : }
2643 : :
2644 [ + + ]: 2764 : if (g_inputdeck.get< tag::pref, tag::pref >())
2645 [ + - ][ + - ]: 402 : elemfieldnames.push_back( "NDOF" );
2646 : :
2647 [ + - ][ + - ]: 2764 : elemfieldnames.push_back( "shock_marker" );
2648 : :
2649 [ - + ][ - - ]: 2764 : Assert( elemfieldnames.size() == m_elemfields.size(), "Size mismatch" );
[ - - ][ - - ]
2650 [ - + ][ - - ]: 2764 : Assert( nodefieldnames.size() == m_nodefields.size(), "Size mismatch" );
[ - - ][ - - ]
2651 : :
2652 : : // Output chare mesh and fields metadata to file
2653 : 2764 : const auto& triinpoel = m_outmesh.triinpoel;
2654 [ + - ]: 8292 : d->write( inpoel, m_outmesh.coord, m_outmesh.bface, {},
2655 [ + - ]: 5528 : tk::remap( triinpoel, lid ), elemfieldnames, nodefieldnames,
2656 : 2764 : {}, m_elemfields, m_nodefields, {}, c );
2657 : 2764 : }
2658 : :
2659 : : void
2660 : 18686 : DG::comnodeout( const std::vector< std::size_t >& gid,
2661 : : const std::vector< std::size_t >& nesup,
2662 : : const std::vector< std::vector< tk::real > >& L )
2663 : : // *****************************************************************************
2664 : : // Receive chare-boundary nodal solution (for field output) contributions from
2665 : : // neighboring chares
2666 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
2667 : : //! \param[in] nesup Number of elements surrounding points
2668 : : //! \param[in] L Partial contributions of node fields to chare-boundary nodes
2669 : : // *****************************************************************************
2670 : : {
2671 [ - + ][ - - ]: 18686 : Assert( gid.size() == nesup.size(), "Size mismatch" );
[ - - ][ - - ]
2672 [ - + ]: 18686 : for (std::size_t f=0; f<L.size(); ++f)
2673 [ - - ][ - - ]: 0 : Assert( gid.size() == L[f].size(), "Size mismatch" );
[ - - ][ - - ]
2674 : :
2675 [ + + ]: 179040 : for (std::size_t i=0; i<gid.size(); ++i) {
2676 : 160354 : auto& nf = m_nodefieldsc[ gid[i] ];
2677 : 160354 : nf.first.resize( L.size() );
2678 [ - + ]: 160354 : for (std::size_t f=0; f<L.size(); ++f) nf.first[f] += L[f][i];
2679 : 160354 : nf.second += nesup[i];
2680 : : }
2681 : :
2682 : : // When we have heard from all chares we communicate with, this chare is done
2683 [ + + ]: 18686 : if (++m_nnod == Disc()->NodeCommMap().size()) {
2684 : 2564 : m_nnod = 0;
2685 : 2564 : comnodeout_complete();
2686 : : }
2687 : 18686 : }
2688 : :
2689 : : void
2690 : 74580 : DG::stage()
2691 : : // *****************************************************************************
2692 : : // Evaluate whether to continue with next time step stage
2693 : : // *****************************************************************************
2694 : : {
2695 : : // Increment Runge-Kutta stage counter
2696 : 74580 : ++m_stage;
2697 : :
2698 : : // if not all Runge-Kutta stages complete, continue to next time stage,
2699 : : // otherwise prepare for nodal field output
2700 [ + + ]: 74580 : if (m_stage < 3)
2701 : 49720 : next();
2702 : : else
2703 [ + - ][ + - ]: 24860 : startFieldOutput( CkCallback(CkIndex_DG::step(), thisProxy[thisIndex]) );
[ + - ]
2704 : 74580 : }
2705 : :
2706 : : void
2707 : 23942 : DG::evalLB( int nrestart )
2708 : : // *****************************************************************************
2709 : : // Evaluate whether to do load balancing
2710 : : //! \param[in] nrestart Number of times restarted
2711 : : // *****************************************************************************
2712 : : {
2713 : 23942 : auto d = Disc();
2714 : :
2715 : : // Detect if just returned from a checkpoint and if so, zero timers
2716 : 23942 : d->restarted( nrestart );
2717 : :
2718 : 23942 : const auto lbfreq = g_inputdeck.get< tag::cmd, tag::lbfreq >();
2719 : 23942 : const auto nonblocking = g_inputdeck.get< tag::cmd, tag::nonblocking >();
2720 : :
2721 : : // Load balancing if user frequency is reached or after the second time-step
2722 [ + + ][ + + ]: 23942 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
[ + + ]
2723 : :
2724 : 20822 : AtSync();
2725 [ - + ]: 20822 : if (nonblocking) next();
2726 : :
2727 : : } else {
2728 : :
2729 : 3120 : next();
2730 : :
2731 : : }
2732 : 23942 : }
2733 : :
2734 : : void
2735 : 23942 : DG::evalRestart()
2736 : : // *****************************************************************************
2737 : : // Evaluate whether to save checkpoint/restart
2738 : : // *****************************************************************************
2739 : : {
2740 : 23942 : auto d = Disc();
2741 : :
2742 : 23942 : const auto rsfreq = g_inputdeck.get< tag::cmd, tag::rsfreq >();
2743 : 23942 : const auto benchmark = g_inputdeck.get< tag::cmd, tag::benchmark >();
2744 : :
2745 [ + + ][ - + ]: 23942 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
[ - + ]
2746 : :
2747 [ - - ]: 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
2748 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
2749 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
2750 : :
2751 : : } else {
2752 : :
2753 : 23942 : evalLB( /* nrestart = */ -1 );
2754 : :
2755 : : }
2756 : 23942 : }
2757 : :
2758 : : void
2759 : 24860 : DG::step()
2760 : : // *****************************************************************************
2761 : : // Evaluate wether to continue with next time step
2762 : : // *****************************************************************************
2763 : : {
2764 : : // Free memory storing output mesh
2765 : 24860 : m_outmesh.destroy();
2766 : :
2767 : 24860 : auto d = Disc();
2768 : :
2769 : : // Output one-liner status report to screen
2770 : 24860 : d->status();
2771 : : // Reset Runge-Kutta stage counter
2772 : 24860 : m_stage = 0;
2773 : :
2774 : 24860 : const auto term = g_inputdeck.get< tag::discr, tag::term >();
2775 : 24860 : const auto nstep = g_inputdeck.get< tag::discr, tag::nstep >();
2776 : 24860 : const auto eps = std::numeric_limits< tk::real >::epsilon();
2777 : :
2778 : : // If neither max iterations nor max time reached, continue, otherwise finish
2779 [ + - ][ + + ]: 24860 : if (std::fabs(d->T()-term) > eps && d->It() < nstep) {
[ + + ]
2780 : :
2781 : 23942 : evalRestart();
2782 : :
2783 : : } else {
2784 : :
2785 : 918 : auto meshid = d->MeshId();
2786 [ + - ]: 918 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
2787 [ + - ][ + - ]: 1836 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
2788 : :
2789 : : }
2790 : 24860 : }
2791 : :
2792 : : #include "NoWarning/dg.def.h"
|
