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 [ + - ]: 5139 : 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 >& interface,
139 : : const std::vector< std::size_t >& ndof );
140 : :
141 : : //! Receive contributions to nodal gradients on chare-boundaries
142 : : void
143 : : comnodalExtrema( const std::vector< std::size_t >& gid,
144 : : const std::vector< std::vector< tk::real > >& G1,
145 : : const std::vector< std::vector< tk::real > >& G2 );
146 : :
147 : : //! Initialize the vector of nodal extrema
148 : : void resizeNodalExtremac();
149 : :
150 : : //! Compute the nodal extrema of ref el derivatives for chare-boundary nodes
151 : : void evalNodalExtrmRefEl(
152 : : const std::size_t ncomp,
153 : : const std::size_t nprim,
154 : : const std::size_t ndof_NodalExtrm,
155 : : const std::vector< std::size_t >& bndel,
156 : : const std::vector< std::size_t >& inpoel,
157 : : const std::vector< std::size_t >& gid,
158 : : const std::unordered_map< std::size_t, std::size_t >& bid,
159 : : const tk::Fields& U,
160 : : const tk::Fields& P,
161 : : std::vector< std::vector<tk::real> >& uNodalExtrm,
162 : : std::vector< std::vector<tk::real> >& pNodalExtrm );
163 : :
164 : : //! \brief Receive nodal solution (ofor field output) contributions from
165 : : //! neighboring chares
166 : : void comnodeout( const std::vector< std::size_t >& gid,
167 : : const std::vector< std::size_t >& nesup,
168 : : const std::vector< std::vector< tk::real > >& Lu,
169 : : const std::vector< std::vector< tk::real > >& Lp );
170 : :
171 : : //! Optionally refine/derefine mesh
172 : : void refine( const std::vector< tk::real >& l2res );
173 : :
174 : : //! Receive new mesh from Refiner
175 : : void resizePostAMR(
176 : : const std::vector< std::size_t >& /* ginpoel */,
177 : : const tk::UnsMesh::Chunk& chunk,
178 : : const tk::UnsMesh::Coords& coord,
179 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /* addedNodes */,
180 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
181 : : const std::set< std::size_t >& removedNodes,
182 : : const std::unordered_map< std::size_t, std::size_t >& amrNodeMap,
183 : : const tk::NodeCommMap& nodeCommMap,
184 : : const std::map< int, std::vector< std::size_t > >& bface,
185 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
186 : : const std::vector< std::size_t >& triinpoel,
187 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
188 : : elemblockid );
189 : :
190 : : //! Extract field output going to file
191 : : void extractFieldOutput(
192 : : const std::vector< std::size_t >& /* ginpoel */,
193 : : const tk::UnsMesh::Chunk& chunk,
194 : : const tk::UnsMesh::Coords& coord,
195 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /* addedNodes */,
196 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
197 : : const tk::NodeCommMap& nodeCommMap,
198 : : const std::map< int, std::vector< std::size_t > >& bface,
199 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
200 : : const std::vector< std::size_t >& triinpoel,
201 : : CkCallback c );
202 : :
203 : : //! Const-ref access to current solution
204 : : //! \return Const-ref to current solution
205 : 0 : const tk::Fields& solution() const { return m_u; }
206 : :
207 : : //! Compute left hand side
208 : : void lhs();
209 : :
210 : : //! Unused in DG
211 : 0 : void resized() {}
212 : :
213 : : //! (no-op)
214 : 0 : void transferSol() {}
215 : :
216 : : //! (no-op)
217 : 0 : void advance( tk::real, tk::real ) {}
218 : :
219 : : //! Compute right hand side and solve system
220 : : void solve( tk::real newdt );
221 : :
222 : : //! Evaluate whether to continue with next time step
223 : : void step();
224 : :
225 : : /** @name Charm++ pack/unpack serializer member functions */
226 : : ///@{
227 : : //! \brief Pack/Unpack serialize member function
228 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
229 : 15926 : void pup( PUP::er &p ) override {
230 : 15926 : p | m_disc;
231 : 15926 : p | m_ghosts;
232 : 15926 : p | m_ndof_NodalExtrm;
233 : 15926 : p | m_nsol;
234 : 15926 : p | m_ninitsol;
235 : 15926 : p | m_nlim;
236 : 15926 : p | m_nnod;
237 : 15926 : p | m_nrefine;
238 : 15926 : p | m_nsmooth;
239 : 15926 : p | m_nreco;
240 : 15926 : p | m_nnodalExtrema;
241 : 15926 : p | m_u;
242 : 15926 : p | m_un;
243 : 15926 : p | m_p;
244 : 15926 : p | m_geoElem;
245 : 15926 : p | m_lhs;
246 : 15926 : p | m_mtInv;
247 : 15926 : p | m_uNodalExtrm;
248 : 15926 : p | m_pNodalExtrm;
249 : 15926 : p | m_uNodalExtrmc;
250 : 15926 : p | m_pNodalExtrmc;
251 : 15926 : p | m_rhs;
252 : 15926 : p | m_rhsprev;
253 : 15926 : p | m_stiffrhs;
254 : 15926 : p | m_stiffrhsprev;
255 : 15926 : p | m_stiffEqIdx;
256 : 15926 : p | m_nonStiffEqIdx;
257 : 15926 : p | m_nstiffeq;
258 : 15926 : p | m_nnonstiffeq;
259 : 15926 : p | m_npoin;
260 : 15926 : p | m_diag;
261 : 15926 : p | m_stage;
262 : 15926 : p | m_ndof;
263 : 15926 : p | m_interface;
264 : 15926 : p | m_numEqDof;
265 : 15926 : p | m_uc;
266 : 15926 : p | m_pc;
267 : 15926 : p | m_ndofc;
268 : 15926 : p | m_interfacec;
269 : 15926 : p | m_initial;
270 : 15926 : p | m_uElemfields;
271 : 15926 : p | m_pElemfields;
272 : 15926 : p | m_uNodefields;
273 : 15926 : p | m_pNodefields;
274 : 15926 : p | m_uNodefieldsc;
275 : 15926 : p | m_pNodefieldsc;
276 : 15926 : p | m_outmesh;
277 : 15926 : p | m_boxelems;
278 : 15926 : p | m_shockmarker;
279 : 15926 : }
280 : : //! \brief Pack/Unpack serialize operator|
281 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
282 : : //! \param[in,out] i DG object reference
283 : : friend void operator|( PUP::er& p, DG& i ) { i.pup(p); }
284 : : //@}
285 : :
286 : : private:
287 : : //! Discretization proxy
288 : : CProxy_Discretization m_disc;
289 : : //! Distributed Ghosts proxy
290 : : CProxy_Ghosts m_ghosts;
291 : : //! \brief Degree of freedom for nodal extrema vector. When DGP1 is applied,
292 : : //! there is one degree of freedom for cell average variable. When DGP2 is
293 : : //! applied, the degree of freedom is 4 which refers to cell average and
294 : : //! gradients in three directions
295 : : std::size_t m_ndof_NodalExtrm;
296 : : //! Counter signaling that we have received all our solution ghost data
297 : : std::size_t m_nsol;
298 : : //! \brief Counter signaling that we have received all our solution ghost
299 : : //! data during setup
300 : : std::size_t m_ninitsol;
301 : : //! \brief Counter signaling that we have received all our limiter function
302 : : //! ghost data
303 : : std::size_t m_nlim;
304 : : //! \brief Counter signaling that we have received all our node solution
305 : : //! contributions
306 : : std::size_t m_nnod;
307 : : //! Counter signaling that we have received all refined ndof
308 : : std::size_t m_nrefine;
309 : : //! Counter signaling that we have received all smoothed ndof
310 : : std::size_t m_nsmooth;
311 : : //! Counter signaling that we have received all our reconstructed ghost data
312 : : std::size_t m_nreco;
313 : : //! \brief Counter signaling that we have received all our nodal extrema from
314 : : //! ghost chare partitions
315 : : std::size_t m_nnodalExtrema;
316 : : //! Counters signaling how many stiff and non-stiff equations in the system
317 : : std::size_t m_nstiffeq, m_nnonstiffeq;
318 : : //! Vector of unknown/solution average over each mesh element
319 : : tk::Fields m_u;
320 : : //! Vector of unknown at previous time-step
321 : : tk::Fields m_un;
322 : : //! Vector of primitive quantities over each mesh element
323 : : tk::Fields m_p;
324 : : //! Element geometry
325 : : tk::Fields m_geoElem;
326 : : //! Left-hand side mass-matrix which is a diagonal matrix
327 : : tk::Fields m_lhs;
328 : : //! Vector of right-hand side
329 : : tk::Fields m_rhs;
330 : : //! Vector of previous right-hand side values used in the IMEX-RK scheme
331 : : tk::Fields m_rhsprev;
332 : : //! Vector of right-hand side for stiff equations
333 : : tk::Fields m_stiffrhs;
334 : : //! Vector of previous right-hand side values for stiff equations
335 : : tk::Fields m_stiffrhsprev;
336 : : //! Vectors that indicates which equations are stiff and non-stiff
337 : : std::vector< std::size_t > m_stiffEqIdx, m_nonStiffEqIdx;
338 : : //! Inverse of Taylor mass-matrix
339 : : std::vector< std::vector< tk::real > > m_mtInv;
340 : : //! Vector of nodal extrema for conservative variables
341 : : std::vector< std::vector<tk::real> > m_uNodalExtrm;
342 : : //! Vector of nodal extrema for primitive variables
343 : : std::vector< std::vector<tk::real> > m_pNodalExtrm;
344 : : //! Buffer for vector of nodal extrema for conservative variables
345 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_uNodalExtrmc;
346 : : //! Buffer for vector of nodal extrema for primitive variables
347 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_pNodalExtrmc;
348 : : //! Counter for number of nodes on this chare excluding ghosts
349 : : std::size_t m_npoin;
350 : : //! Diagnostics object
351 : : ElemDiagnostics m_diag;
352 : : //! Runge-Kutta stage counter
353 : : std::size_t m_stage;
354 : : //! Vector of local number of degrees of freedom for each element
355 : : std::vector< std::size_t > m_ndof;
356 : : //! Interface marker for field output
357 : : std::vector< std::size_t > m_interface;
358 : : //! Vector of number of degrees of freedom for each PDE equation/component
359 : : std::vector< std::size_t > m_numEqDof;
360 : : //! Solution receive buffers for ghosts only
361 : : std::array< std::vector< std::vector< tk::real > >, 3 > m_uc;
362 : : //! Primitive-variable receive buffers for ghosts only
363 : : std::array< std::vector< std::vector< tk::real > >, 3 > m_pc;
364 : : //! \brief Number of degrees of freedom (for p-adaptive) receive buffers
365 : : //! for ghosts only
366 : : std::array< std::vector< std::size_t >, 3 > m_ndofc;
367 : : //! Interface marker receive buffers for ghost only
368 : : std::array< std::vector< std::size_t >, 1 > m_interfacec;
369 : : //! 1 if starting time stepping, 0 if during time stepping
370 : : std::size_t m_initial;
371 : : //! Solution elem output fields
372 : : tk::Fields m_uElemfields;
373 : : //! Primitive elem output fields
374 : : tk::Fields m_pElemfields;
375 : : //! Solution nodal output fields
376 : : tk::Fields m_uNodefields;
377 : : //! Primitive nodal output fields
378 : : tk::Fields m_pNodefields;
379 : : //! Receive buffer for communication of solution node fields
380 : : //! \details Key: global node id, value: output fields and number of
381 : : //! elements surrounding the node
382 : : std::unordered_map< std::size_t, std::pair< std::vector< tk::real >,
383 : : std::size_t > > m_uNodefieldsc;
384 : : //! Receive buffer for communication of primitive quantity node fields
385 : : //! \details Key: global node id, value: output fields and number of
386 : : //! elements surrounding the node
387 : : std::unordered_map< std::size_t, std::pair< std::vector< tk::real >,
388 : : std::size_t > > m_pNodefieldsc;
389 : : //! Storage for refined mesh used for field output
390 : : Ghosts::OutMesh m_outmesh;
391 : : //! Element ids at which box ICs are defined by user (multiple boxes)
392 : : std::vector< std::unordered_set< std::size_t > > m_boxelems;
393 : : //! Shock detection marker for field output
394 : : std::vector< std::size_t > m_shockmarker;
395 : :
396 : : //! Access bound Discretization class pointer
397 : 1090937 : Discretization* Disc() const {
398 [ + - ][ - + ]: 1090937 : Assert( m_disc[ thisIndex ].ckLocal() != nullptr, "ckLocal() null" );
[ - - ][ - - ]
[ - - ]
399 [ + - ]: 1090937 : return m_disc[ thisIndex ].ckLocal();
400 : : }
401 : :
402 : : //! Access bound Discretization class pointer
403 : 188291699 : Ghosts* myGhosts() const {
404 [ + - ][ - + ]: 188291699 : Assert( m_ghosts[ thisIndex ].ckLocal() != nullptr, "ckLocal() null" );
[ - - ][ - - ]
[ - - ]
405 [ + - ]: 188291699 : return m_ghosts[ thisIndex ].ckLocal();
406 : : }
407 : :
408 : : //! Output mesh field data
409 : : void writeFields(
410 : : CkCallback c,
411 : : const std::unordered_map< std::size_t, std::size_t >& addedTets );
412 : :
413 : : //! Add the protective layer for ndof refinement
414 : : void refine();
415 : :
416 : : //! Smooth the refined ndof distribution
417 : : void smooth();
418 : :
419 : : //! Compute solution reconstructions
420 : : void reco();
421 : :
422 : : //! Compute nodal extrema at chare-boundary nodes
423 : : void nodalExtrema();
424 : :
425 : : //! Compute limiter function
426 : : void lim();
427 : :
428 : : //! Compute time step size
429 : : void dt();
430 : :
431 : : //! Evaluate whether to continue with next time step stage
432 : : void stage();
433 : :
434 : : //! Evaluate whether to save checkpoint/restart
435 : : void evalRestart();
436 : :
437 : : //! p-refine all elements that are adjacent to p-refined elements
438 : : void refine_ndof();
439 : :
440 : : //! Smooth the refined ndof distribution to avoid zigzag refinement
441 : : void smooth_ndof();
442 : :
443 : : //! Decide wether to output field data
444 : : bool fieldOutput() const;
445 : :
446 : : //! Decide if we write field output using a refined mesh
447 : : bool refinedOutput() const;
448 : :
449 : : //! Start preparing fields for output to file
450 : : void startFieldOutput( CkCallback c );
451 : :
452 : : //! Compute the integration step for IMEX-RK
453 : : void imex_integrate();
454 : : };
455 : :
456 : : } // inciter::
457 : :
458 : : #endif // DG_h
|
