Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/DG.hpp
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 : :
13 : : There are a potentially large number of DG Charm++ chares created by
14 : : Transporter. Each DG gets a chunk of the full load (part of the mesh)
15 : : and does the same: initializes and advances a number of PDE systems in time.
16 : :
17 : : The implementation uses the Charm++ runtime system and is fully
18 : : asynchronous, overlapping computation and communication. The algorithm
19 : : utilizes the structured dagger (SDAG) Charm++ functionality. The high-level
20 : : overview of the algorithm structure and how it interfaces with Charm++ is
21 : : discussed in the Charm++ interface file src/Inciter/dg.ci.
22 : : */
23 : : // *****************************************************************************
24 : : #ifndef DG_h
25 : : #define DG_h
26 : :
27 : : #include <array>
28 : : #include <set>
29 : : #include <unordered_set>
30 : : #include <unordered_map>
31 : :
32 : : #include "DerivedData.hpp"
33 : : #include "FaceData.hpp"
34 : : #include "ElemDiagnostics.hpp"
35 : : #include "Ghosts.hpp"
36 : :
37 : : #include "NoWarning/dg.decl.h"
38 : :
39 : : namespace inciter {
40 : :
41 : : //! DG Charm++ chare array used to advance PDEs in time with DG+RK
42 : : class DG : public CBase_DG {
43 : :
44 : : public:
45 : : #if defined(__clang__)
46 : : #pragma clang diagnostic push
47 : : #pragma clang diagnostic ignored "-Wunused-parameter"
48 : : #pragma clang diagnostic ignored "-Wdeprecated-declarations"
49 : : #elif defined(STRICT_GNUC)
50 : : #pragma GCC diagnostic push
51 : : #pragma GCC diagnostic ignored "-Wunused-parameter"
52 : : #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
53 : : #elif defined(__INTEL_COMPILER)
54 : : #pragma warning( push )
55 : : #pragma warning( disable: 1478 )
56 : : #endif
57 : : // Include Charm++ SDAG code. See http://charm.cs.illinois.edu/manuals/html/
58 : : // charm++/manual.html, Sec. "Structured Control Flow: Structured Dagger".
59 : : DG_SDAG_CODE
60 : : #if defined(__clang__)
61 : : #pragma clang diagnostic pop
62 : : #elif defined(STRICT_GNUC)
63 : : #pragma GCC diagnostic pop
64 : : #elif defined(__INTEL_COMPILER)
65 : : #pragma warning( pop )
66 : : #endif
67 : :
68 : : //! Constructor
69 : : explicit DG( const CProxy_Discretization& disc,
70 : : const CProxy_Ghosts& ghostsproxy,
71 : : const std::map< int, std::vector< std::size_t > >& bface,
72 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
73 : : const std::vector< std::size_t >& triinpoel );
74 : :
75 : : #if defined(__clang__)
76 : : #pragma clang diagnostic push
77 : : #pragma clang diagnostic ignored "-Wundefined-func-template"
78 : : #endif
79 : : //! Migrate constructor
80 : : // cppcheck-suppress uninitMemberVar
81 : 10182 : explicit DG( CkMigrateMessage* msg ) : CBase_DG( msg ) {}
82 : : #if defined(__clang__)
83 : : #pragma clang diagnostic pop
84 : : #endif
85 : :
86 : : //! Return from migration
87 : : void ResumeFromSync() override;
88 : :
89 : : //! Configure Charm++ reduction types for concatenating BC nodelists
90 : : static void registerReducers();
91 : :
92 : : //! Resize solution vectors after extension due to Ghosts and continue setup
93 : : void resizeSolVectors();
94 : :
95 : : //! Setup: query boundary conditions, output mesh, etc.
96 : : void setup();
97 : :
98 : : //! Receive total box IC volume and set conditions in box
99 : : void box( tk::real v, const std::vector< tk::real >& blkvols );
100 : :
101 : : // Evaluate whether to do load balancing
102 : : void evalLB( int nrestart );
103 : :
104 : : //! Start time stepping
105 : : void start();
106 : :
107 : : //! Continue to next time step
108 : : void next();
109 : :
110 : : //! Receive chare-boundary updated solution from neighboring chares
111 : : void comrefine( int fromch,
112 : : const std::vector< std::size_t >& tetid,
113 : : const std::vector< std::size_t >& ndof );
114 : :
115 : : // Receive chare-boundary refined ndof from neighboring chares
116 : : void comsmooth( int fromch,
117 : : const std::vector< std::size_t >& tetid,
118 : : const std::vector< std::size_t >& ndof );
119 : :
120 : : //! Receive chare-boundary limiter function data from neighboring chares
121 : : void comlim( int fromch,
122 : : const std::vector< std::size_t >& tetid,
123 : : const std::vector< std::vector< tk::real > >& u,
124 : : const std::vector< std::vector< tk::real > >& prim );
125 : :
126 : : //! Receive chare-boundary reconstructed data from neighboring chares
127 : : void comreco( int fromch,
128 : : const std::vector< std::size_t >& tetid,
129 : : const std::vector< std::vector< tk::real > >& u,
130 : : const std::vector< std::vector< tk::real > >& prim );
131 : :
132 : : //! Receive chare-boundary ghost data from neighboring chares
133 : : void comsol( int fromch,
134 : : std::size_t fromstage,
135 : : const std::vector< std::size_t >& tetid,
136 : : const std::vector< std::vector< tk::real > >& u,
137 : : const std::vector< std::vector< tk::real > >& prim,
138 : : const std::vector< std::size_t >& ndof );
139 : :
140 : : //! Receive contributions to nodal gradients on chare-boundaries
141 : : void
142 : : comnodalExtrema( const std::vector< std::size_t >& gid,
143 : : const std::vector< std::vector< tk::real > >& G1,
144 : : const std::vector< std::vector< tk::real > >& G2 );
145 : :
146 : : //! Initialize the vector of nodal extrema
147 : : void resizeNodalExtremac();
148 : :
149 : : //! Compute the nodal extrema of ref el derivatives for chare-boundary nodes
150 : : void evalNodalExtrmRefEl(
151 : : const std::size_t ncomp,
152 : : const std::size_t nprim,
153 : : const std::size_t ndof_NodalExtrm,
154 : : const std::vector< std::size_t >& bndel,
155 : : const std::vector< std::size_t >& inpoel,
156 : : const std::vector< std::size_t >& gid,
157 : : const std::unordered_map< std::size_t, std::size_t >& bid,
158 : : const tk::Fields& U,
159 : : const tk::Fields& P,
160 : : std::vector< std::vector<tk::real> >& uNodalExtrm,
161 : : std::vector< std::vector<tk::real> >& pNodalExtrm );
162 : :
163 : : //! \brief Receive nodal solution (ofor field output) contributions from
164 : : //! neighboring chares
165 : : void comnodeout( const std::vector< std::size_t >& gid,
166 : : const std::vector< std::size_t >& nesup,
167 : : const std::vector< std::vector< tk::real > >& Lu,
168 : : const std::vector< std::vector< tk::real > >& Lp );
169 : :
170 : : //! Optionally refine/derefine mesh
171 : : void refine( const std::vector< tk::real >& l2res );
172 : :
173 : : //! Receive new mesh from Refiner
174 : : void resizePostAMR(
175 : : const std::vector< std::size_t >& /* ginpoel */,
176 : : const tk::UnsMesh::Chunk& chunk,
177 : : const tk::UnsMesh::Coords& coord,
178 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /* addedNodes */,
179 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
180 : : const std::set< std::size_t >& removedNodes,
181 : : const std::unordered_map< std::size_t, std::size_t >& amrNodeMap,
182 : : const tk::NodeCommMap& nodeCommMap,
183 : : const std::map< int, std::vector< std::size_t > >& bface,
184 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
185 : : const std::vector< std::size_t >& triinpoel,
186 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
187 : : elemblockid );
188 : :
189 : : //! Extract field output going to file
190 : : void extractFieldOutput(
191 : : const std::vector< std::size_t >& /* ginpoel */,
192 : : const tk::UnsMesh::Chunk& chunk,
193 : : const tk::UnsMesh::Coords& coord,
194 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /* addedNodes */,
195 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
196 : : const tk::NodeCommMap& nodeCommMap,
197 : : const std::map< int, std::vector< std::size_t > >& bface,
198 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
199 : : const std::vector< std::size_t >& triinpoel,
200 : : CkCallback c );
201 : :
202 : : //! Const-ref access to current solution
203 : : //! \return Const-ref to current solution
204 : : const tk::Fields& solution() const { return m_u; }
205 : :
206 : : //! Compute left hand side
207 : : void lhs();
208 : :
209 : : //! Unused in DG
210 : : void resized() {}
211 : :
212 : : //! (no-op)
213 : : void transferSol() {}
214 : :
215 : : //! (no-op)
216 : : void advance( tk::real, tk::real ) {}
217 : :
218 : : //! Compute right hand side and solve system
219 : : void solve( tk::real newdt );
220 : :
221 : : //! Evaluate whether to continue with next time step
222 : : void step();
223 : :
224 : : /** @name Charm++ pack/unpack serializer member functions */
225 : : ///@{
226 : : //! \brief Pack/Unpack serialize member function
227 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
228 : 15834 : void pup( PUP::er &p ) override {
229 : 15834 : p | m_disc;
230 : 15834 : p | m_ghosts;
231 : 15834 : p | m_ndof_NodalExtrm;
232 : 15834 : p | m_nsol;
233 : 15834 : p | m_ninitsol;
234 : 15834 : p | m_nlim;
235 : 15834 : p | m_nnod;
236 : 15834 : p | m_nrefine;
237 : 15834 : p | m_nsmooth;
238 : 15834 : p | m_nreco;
239 : 15834 : p | m_nnodalExtrema;
240 : 15834 : p | m_u;
241 : 15834 : p | m_un;
242 : 15834 : p | m_p;
243 : 15834 : p | m_geoElem;
244 : 15834 : p | m_lhs;
245 : 15834 : p | m_mtInv;
246 : 15834 : p | m_uNodalExtrm;
247 : 15834 : p | m_pNodalExtrm;
248 : 15834 : p | m_uNodalExtrmc;
249 : 15834 : p | m_pNodalExtrmc;
250 : 15834 : p | m_rhs;
251 : 15834 : p | m_npoin;
252 : : p | m_diag;
253 : 15834 : p | m_stage;
254 : 15834 : p | m_ndof;
255 : 15834 : p | m_numEqDof;
256 : : p | m_uc;
257 : : p | m_pc;
258 : : p | m_ndofc;
259 : : p | m_initial;
260 : 15834 : p | m_uElemfields;
261 : 15834 : p | m_pElemfields;
262 : 15834 : p | m_uNodefields;
263 : 15834 : p | m_pNodefields;
264 : 15834 : p | m_uNodefieldsc;
265 : 15834 : p | m_pNodefieldsc;
266 : 15834 : p | m_outmesh;
267 : 15834 : p | m_boxelems;
268 : 15834 : p | m_shockmarker;
269 : 15834 : }
270 : : //! \brief Pack/Unpack serialize operator|
271 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
272 : : //! \param[in,out] i DG object reference
273 : : friend void operator|( PUP::er& p, DG& i ) { i.pup(p); }
274 : : //@}
275 : :
276 : : private:
277 : : //! Discretization proxy
278 : : CProxy_Discretization m_disc;
279 : : //! Distributed Ghosts proxy
280 : : CProxy_Ghosts m_ghosts;
281 : : //! \brief Degree of freedom for nodal extrema vector. When DGP1 is applied,
282 : : //! there is one degree of freedom for cell average variable. When DGP2 is
283 : : //! applied, the degree of freedom is 4 which refers to cell average and
284 : : //! gradients in three directions
285 : : std::size_t m_ndof_NodalExtrm;
286 : : //! Counter signaling that we have received all our solution ghost data
287 : : std::size_t m_nsol;
288 : : //! \brief Counter signaling that we have received all our solution ghost
289 : : //! data during setup
290 : : std::size_t m_ninitsol;
291 : : //! \brief Counter signaling that we have received all our limiter function
292 : : //! ghost data
293 : : std::size_t m_nlim;
294 : : //! \brief Counter signaling that we have received all our node solution
295 : : //! contributions
296 : : std::size_t m_nnod;
297 : : //! Counter signaling that we have received all refined ndof
298 : : std::size_t m_nrefine;
299 : : //! Counter signaling that we have received all smoothed ndof
300 : : std::size_t m_nsmooth;
301 : : //! Counter signaling that we have received all our reconstructed ghost data
302 : : std::size_t m_nreco;
303 : : //! \brief Counter signaling that we have received all our nodal extrema from
304 : : //! ghost chare partitions
305 : : std::size_t m_nnodalExtrema;
306 : : //! Vector of unknown/solution average over each mesh element
307 : : tk::Fields m_u;
308 : : //! Vector of unknown at previous time-step
309 : : tk::Fields m_un;
310 : : //! Vector of primitive quantities over each mesh element
311 : : tk::Fields m_p;
312 : : //! Element geometry
313 : : tk::Fields m_geoElem;
314 : : //! Left-hand side mass-matrix which is a diagonal matrix
315 : : tk::Fields m_lhs;
316 : : //! Vector of right-hand side
317 : : tk::Fields m_rhs;
318 : : //! Inverse of Taylor mass-matrix
319 : : std::vector< std::vector< tk::real > > m_mtInv;
320 : : //! Vector of nodal extrema for conservative variables
321 : : std::vector< std::vector<tk::real> > m_uNodalExtrm;
322 : : //! Vector of nodal extrema for primitive variables
323 : : std::vector< std::vector<tk::real> > m_pNodalExtrm;
324 : : //! Buffer for vector of nodal extrema for conservative variables
325 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_uNodalExtrmc;
326 : : //! Buffer for vector of nodal extrema for primitive variables
327 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_pNodalExtrmc;
328 : : //! Counter for number of nodes on this chare excluding ghosts
329 : : std::size_t m_npoin;
330 : : //! Diagnostics object
331 : : ElemDiagnostics m_diag;
332 : : //! Runge-Kutta stage counter
333 : : std::size_t m_stage;
334 : : //! Vector of local number of degrees of freedom for each element
335 : : std::vector< std::size_t > m_ndof;
336 : : //! Vector of number of degrees of freedom for each PDE equation/component
337 : : std::vector< std::size_t > m_numEqDof;
338 : : //! Solution receive buffers for ghosts only
339 : : std::array< std::vector< std::vector< tk::real > >, 3 > m_uc;
340 : : //! Primitive-variable receive buffers for ghosts only
341 : : std::array< std::vector< std::vector< tk::real > >, 3 > m_pc;
342 : : //! \brief Number of degrees of freedom (for p-adaptive) receive buffers
343 : : //! for ghosts only
344 : : std::array< std::vector< std::size_t >, 3 > m_ndofc;
345 : : //! 1 if starting time stepping, 0 if during time stepping
346 : : std::size_t m_initial;
347 : : //! Solution elem output fields
348 : : tk::Fields m_uElemfields;
349 : : //! Primitive elem output fields
350 : : tk::Fields m_pElemfields;
351 : : //! Solution nodal output fields
352 : : tk::Fields m_uNodefields;
353 : : //! Primitive nodal output fields
354 : : tk::Fields m_pNodefields;
355 : : //! Receive buffer for communication of solution node fields
356 : : //! \details Key: global node id, value: output fields and number of
357 : : //! elements surrounding the node
358 : : std::unordered_map< std::size_t, std::pair< std::vector< tk::real >,
359 : : std::size_t > > m_uNodefieldsc;
360 : : //! Receive buffer for communication of primitive quantity node fields
361 : : //! \details Key: global node id, value: output fields and number of
362 : : //! elements surrounding the node
363 : : std::unordered_map< std::size_t, std::pair< std::vector< tk::real >,
364 : : std::size_t > > m_pNodefieldsc;
365 : : //! Storage for refined mesh used for field output
366 : : Ghosts::OutMesh m_outmesh;
367 : : //! Element ids at which box ICs are defined by user (multiple boxes)
368 : : std::vector< std::unordered_set< std::size_t > > m_boxelems;
369 : : //! Shock detection marker for field output
370 : : std::vector< std::size_t > m_shockmarker;
371 : :
372 : : //! Access bound Discretization class pointer
373 : 1188642 : Discretization* Disc() const {
374 : : Assert( m_disc[ thisIndex ].ckLocal() != nullptr, "ckLocal() null" );
375 : 2377284 : return m_disc[ thisIndex ].ckLocal();
376 : : }
377 : :
378 : : //! Access bound Discretization class pointer
379 : 84596934 : Ghosts* myGhosts() const {
380 : : Assert( m_ghosts[ thisIndex ].ckLocal() != nullptr, "ckLocal() null" );
381 : 169193868 : return m_ghosts[ thisIndex ].ckLocal();
382 : : }
383 : :
384 : : //! Output mesh field data
385 : : void writeFields(
386 : : CkCallback c,
387 : : const std::unordered_map< std::size_t, std::size_t >& addedTets );
388 : :
389 : : //! Add the protective layer for ndof refinement
390 : : void refine();
391 : :
392 : : //! Smooth the refined ndof distribution
393 : : void smooth();
394 : :
395 : : //! Compute solution reconstructions
396 : : void reco();
397 : :
398 : : //! Compute nodal extrema at chare-boundary nodes
399 : : void nodalExtrema();
400 : :
401 : : //! Compute limiter function
402 : : void lim();
403 : :
404 : : //! Compute time step size
405 : : void dt();
406 : :
407 : : //! Evaluate whether to continue with next time step stage
408 : : void stage();
409 : :
410 : : //! Evaluate whether to save checkpoint/restart
411 : : void evalRestart();
412 : :
413 : : //! p-refine all elements that are adjacent to p-refined elements
414 : : void refine_ndof();
415 : :
416 : : //! Smooth the refined ndof distribution to avoid zigzag refinement
417 : : void smooth_ndof();
418 : :
419 : : //! Decide wether to output field data
420 : : bool fieldOutput() const;
421 : :
422 : : //! Decide if we write field output using a refined mesh
423 : : bool refinedOutput() const;
424 : :
425 : : //! Start preparing fields for output to file
426 : : void startFieldOutput( CkCallback c );
427 : : };
428 : :
429 : : } // inciter::
430 : :
431 : : #endif // DG_h
|
