Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/ZalCG.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.,
7 : : 2022-2024 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief ZalCG: Taylor-Galerkin, FCT, edge-based continuous Galerkin
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #pragma once
14 : :
15 : : #include <vector>
16 : : #include <map>
17 : :
18 : : #include "Types.hpp"
19 : : #include "Fields.hpp"
20 : : #include "Table.hpp"
21 : : #include "DerivedData.hpp"
22 : : #include "NodeDiagnostics.hpp"
23 : : #include "InciterConfig.hpp"
24 : :
25 : : #include "NoWarning/zalcg.decl.h"
26 : :
27 : : namespace inciter {
28 : :
29 : : extern ctr::Config g_cfg;
30 : :
31 : : //! ZalCG Charm++ chare array used to advance PDEs in time with ZalCG
32 : : class ZalCG : public CBase_ZalCG {
33 : :
34 : : public:
35 : : #if defined(__clang__)
36 : : #pragma clang diagnostic push
37 : : #pragma clang diagnostic ignored "-Wunused-parameter"
38 : : #pragma clang diagnostic ignored "-Wdeprecated-declarations"
39 : : #elif defined(STRICT_GNUC)
40 : : #pragma GCC diagnostic push
41 : : #pragma GCC diagnostic ignored "-Wunused-parameter"
42 : : #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
43 : : #elif defined(__INTEL_COMPILER)
44 : : #pragma warning( push )
45 : : #pragma warning( disable: 1478 )
46 : : #endif
47 : : // Include Charm++ SDAG code. See http://charm.cs.illinois.edu/manuals/html/
48 : : // charm++/manual.html, Sec. "Structured Control Flow: Structured Dagger".
49 : : ZalCG_SDAG_CODE
50 : : #if defined(__clang__)
51 : : #pragma clang diagnostic pop
52 : : #elif defined(STRICT_GNUC)
53 : : #pragma GCC diagnostic pop
54 : : #elif defined(__INTEL_COMPILER)
55 : : #pragma warning( pop )
56 : : #endif
57 : :
58 : : //! Constructor
59 : : explicit ZalCG( const CProxy_Discretization& disc,
60 : : const std::map< int, std::vector< std::size_t > >& bface,
61 : : const std::map< int, std::vector< std::size_t > >& bnode,
62 : : const std::vector< std::size_t >& triinpoel );
63 : :
64 : : #if defined(__clang__)
65 : : #pragma clang diagnostic push
66 : : #pragma clang diagnostic ignored "-Wundefined-func-template"
67 : : #endif
68 : : //! Migrate constructor
69 : : // cppcheck-suppress uninitMemberVar
70 : 897 : explicit ZalCG( CkMigrateMessage* m ) : CBase_ZalCG( m ) {}
71 : : #if defined(__clang__)
72 : : #pragma clang diagnostic pop
73 : : #endif
74 : :
75 : : //! Configure Charm++ custom reduction types initiated from this chare array
76 : : static void registerReducers();
77 : :
78 : : //! Return from migration
79 : : void ResumeFromSync() override;
80 : :
81 : : //! Start setup for solution
82 : : void setup( tk::real v );
83 : :
84 : : // Start time stepping
85 : : void start();
86 : :
87 : : //! Advance equations to next time step
88 : : void advance( tk::real newdt );
89 : :
90 : : //! Start (re-)computing domain and boundary integrals
91 : : void feop();
92 : :
93 : : //! Receive contributions to boundary point normals on chare-boundaries
94 : : void comnorm( const std::unordered_map< int,
95 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > >& inbnd );
96 : :
97 : : //! Receive contributions to right-hand side vector on chare-boundaries
98 : : void comrhs( const std::unordered_map< std::size_t,
99 : : std::vector< tk::real > >& inrhs );
100 : :
101 : : //! Receive antidiffusive and low-order contributions on chare-boundaries
102 : : void comaec( const std::unordered_map< std::size_t,
103 : : std::vector< tk::real > >& inaec );
104 : :
105 : :
106 : : //! Receive allowed limits contributions on chare-boundaries
107 : : void comalw( const std::unordered_map< std::size_t,
108 : : std::vector< tk::real > >& inalw );
109 : :
110 : : //! Receive limited antidiffusive contributions on chare-boundaries
111 : : void comlim( const std::unordered_map< std::size_t,
112 : : std::vector< tk::real > >& inlim );
113 : :
114 : : //! Evaluate residuals
115 : : void evalres( const std::vector< tk::real >& l2res );
116 : :
117 : : //! Communicate deactivation requests
118 : : void comdea( std::size_t reactivate );
119 : :
120 : : //! Receive aggregate deactivation status
121 : : void deastat( int dea );
122 : :
123 : : //! Communicate deactivation status
124 : : void comact( int ch, int deactivated );
125 : :
126 : : //! Receive new mesh from Refiner
127 : : void resizePostAMR(
128 : : const std::vector< std::size_t >& ginpoel,
129 : : const tk::UnsMesh::Chunk& chunk,
130 : : const tk::UnsMesh::Coords& coord,
131 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
132 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
133 : : const std::set< std::size_t >& removedNodes,
134 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
135 : : nodeCommMap,
136 : : const std::map< int, std::vector< std::size_t > >& bface,
137 : : const std::map< int, std::vector< std::size_t > >& bnode,
138 : : const std::vector< std::size_t >& triinpoel );
139 : :
140 : : //! Const-ref access to current solution
141 : : //! \return Const-ref to current solution
142 : : const tk::Fields& solution() const { return m_u; }
143 : :
144 : : //! Compute integral quantities for output
145 : : void integrals();
146 : :
147 : : //! Evaluate whether to continue with next time step
148 : : void step();
149 : :
150 : : // Evaluate whether to do load balancing
151 : : void evalLB( int nrestart );
152 : :
153 : : /** @name Charm++ pack/unpack serializer member functions */
154 : : ///@{
155 : : //! \brief Pack/Unpack serialize member function
156 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
157 : 2829 : void pup( PUP::er &p ) override {
158 : : p | m_disc;
159 : 2829 : p | m_nrhs;
160 : 2829 : p | m_nnorm;
161 : 2829 : p | m_naec;
162 : 2829 : p | m_nalw;
163 : 2829 : p | m_nlim;
164 : 2829 : p | m_ndea;
165 : 2829 : p | m_nact;
166 : 2829 : p | m_todeactivate;
167 : 2829 : p | m_toreactivate;
168 : 2829 : p | m_deactivated;
169 : 2829 : p | m_inactive;
170 : 2829 : p | m_bnode;
171 : 2829 : p | m_bface;
172 : 2829 : p | m_triinpoel;
173 : 2829 : p | m_u;
174 : 2829 : p | m_p;
175 : 2829 : p | m_pc;
176 : 2829 : p | m_q;
177 : 2829 : p | m_qc;
178 : 2829 : p | m_a;
179 : 2829 : p | m_ac;
180 : : // do not pup these, will recompute after migration anyway
181 [ + + ]: 2829 : if (p.isUnpacking()) {
182 : : m_rhs.resize( m_u.nunk(), m_u.nprop() );
183 : : }
184 : 2829 : p | m_rhsc;
185 : 2829 : p | m_vol;
186 : : p | m_diag;
187 : 2829 : p | m_bnorm;
188 : 2829 : p | m_bnormc;
189 : 2829 : p | m_bndpoinint;
190 : 2829 : p | m_bndedgeint;
191 : 2829 : p | m_domedgeint;
192 : 2829 : p | m_bpoin;
193 : 2829 : p | m_bpint;
194 : : p | m_dsupedge;
195 : : p | m_dsupint;
196 : : p | m_dsuplim;
197 : 2829 : p | m_chbndedge;
198 : 2829 : p | m_besym;
199 : 2829 : p | m_dirbcmasks;
200 : 2829 : p | m_prebcnodes;
201 : 2829 : p | m_prebcvals;
202 : 2829 : p | m_symbcnodeset;
203 : 2829 : p | m_symbcnodes;
204 : 2829 : p | m_symbcnorms;
205 : 2829 : p | m_farbcnodeset;
206 : 2829 : p | m_farbcnodes;
207 : 2829 : p | m_farbcnorms;
208 : 2829 : p | m_surfint;
209 : 2829 : p | m_dtp;
210 : 2829 : p | m_tp;
211 : 2829 : p | m_finished;
212 : 2829 : p | m_freezeflow;
213 : 2829 : p | m_fctfreeze;
214 : 2829 : }
215 : : //! \brief Pack/Unpack serialize operator|
216 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
217 : : //! \param[in,out] i ZalCG object reference
218 : : friend void operator|( PUP::er& p, ZalCG& i ) { i.pup(p); }
219 : : //@}
220 : :
221 : : private:
222 : : //! Discretization proxy
223 : : CProxy_Discretization m_disc;
224 : : //! Counter for right-hand side vector nodes updated
225 : : std::size_t m_nrhs;
226 : : //! Counter for receiving boundary point normals
227 : : std::size_t m_nnorm;
228 : : //! Counter for receiving antidiffusive contributions
229 : : std::size_t m_naec;
230 : : //! Counter for receiving allowed limits
231 : : std::size_t m_nalw;
232 : : //! Counter for receiving limited antidiffusive contributions
233 : : std::size_t m_nlim;
234 : : //! Counter for receiving deactivation requests
235 : : std::size_t m_ndea;
236 : : //! Counter for receiving activation status communications
237 : : std::size_t m_nact;
238 : : //! Flag: 1 if chare desires to deactivate
239 : : int m_todeactivate;
240 : : //! Flag: 1 if chare desires to reactivate
241 : : int m_toreactivate;
242 : : //! Flag: 1 if chare is deactivated, 0 if active
243 : : int m_deactivated;
244 : : //! Deactived chares this chare communicates with
245 : : std::unordered_set< int > m_inactive;
246 : : //! Boundary node lists mapped to side set ids used in the input file
247 : : std::map< int, std::vector< std::size_t > > m_bnode;
248 : : //! Boundary face lists mapped to side set ids used in the input file
249 : : std::map< int, std::vector< std::size_t > > m_bface;
250 : : //! Chare-boundary triangle face connecitivity
251 : : std::vector< std::size_t > m_triinpoel;
252 : : //! Unknown/solution vector at mesh nodes
253 : : tk::Fields m_u;
254 : : //! Max/min antidiffusive edge contributions at mesh nodes
255 : : tk::Fields m_p;
256 : : //! Receive buffer for max/min antidiffusive edge contributions
257 : : //! \details Key: global node id, value: max/min antidiff edge contributions
258 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_pc;
259 : : //! Max/min allowed limits at mesh nodes
260 : : tk::Fields m_q;
261 : : //! Receive buffer for max/min allowed limits
262 : : //! \details Key: global node id, value: max/min allowed limits in nodes.
263 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_qc;
264 : : //! Limited antidiffusive contributions at mesh nodes
265 : : tk::Fields m_a;
266 : : //! Receive buffer for limited antidiffusive contributions
267 : : //! \details Key: global node id, value: limited antidiffusive contributions
268 : : //! in nodes.
269 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_ac;
270 : : //! Right-hand side vector (for the high order system)
271 : : tk::Fields m_rhs;
272 : : //! Receive buffer for communication of the right hand side
273 : : //! \details Key: global node id, value: rhs for all scalar components per
274 : : //! node.
275 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_rhsc;
276 : : //! Diagnostics object
277 : : NodeDiagnostics m_diag;
278 : : //! Boundary point normals
279 : : //! \details Outer key: side set id. Inner key: global node id of boundary
280 : : //! point, value: weighted normals, inverse distance square, nodal area.
281 : : std::unordered_map< int,
282 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > > m_bnorm;
283 : : //! Boundary point normals receive buffer
284 : : //! \details Outer key: side set id. Inner key: global node id of boundary
285 : : //! point, value: weighted normals and inverse distance square.
286 : : decltype(m_bnorm) m_bnormc;
287 : : //! Boundary point integrals
288 : : //! \details Key: global node id of boundary point, value: boundary point
289 : : //! integral contributions.
290 : : std::unordered_map< std::size_t, std::array< tk::real, 3 > > m_bndpoinint;
291 : : //! Boundary edge integrals
292 : : //! \details Key: boundary edge-end points with global node ids, value:
293 : : //! boundary edge integral contributions.
294 : : std::unordered_map< tk::UnsMesh::Edge, std::array< tk::real, 3 >,
295 : : tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> > m_bndedgeint;
296 : : //! Domain edge integrals
297 : : std::unordered_map< tk::UnsMesh::Edge, std::array< tk::real, 4 >,
298 : : tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> > m_domedgeint;
299 : : //! Streamable boundary point local ids
300 : : std::vector< std::size_t > m_bpoin;
301 : : //! Streamable boundary point integrals
302 : : std::vector< tk::real > m_bpint;
303 : : //! Superedge (tet, face, edge) end points with local ids for domain edges
304 : : std::array< std::vector< std::size_t >, 3 > m_dsupedge;
305 : : //! Superedge (tet, face, edge) domain edge integrals
306 : : std::array< std::vector< tk::real >, 3 > m_dsupint;
307 : : //! FCT limiter coefficients in domain superedges
308 : : std::array< std::vector< tk::real >, 3 > m_dsuplim;
309 : : //! Chare-boundary edge end-points with difffusion integral
310 : : //! \details Key: neighbor chare id, value: domain-edge end-points and
311 : : //! diffusion integral associated to the edge
312 : : std::unordered_map< int,
313 : : std::vector< std::tuple< tk::UnsMesh::Edge, tk::real > > > m_chbndedge;
314 : : //! Streamable boundary point symmetry BC flags
315 : : std::vector< std::uint8_t > m_besym;
316 : : //! Nodal volumes dynamically adjusted for deactivated chares
317 : : std::vector< tk::real > m_vol;
318 : : //! Nodes and their Dirichlet BC masks
319 : : std::vector< std::size_t > m_dirbcmasks;
320 : : //! Nodes at pressure BCs
321 : : std::vector< std::size_t > m_prebcnodes;
322 : : //! Density and pressure values at pressure BCs
323 : : std::vector< tk::real > m_prebcvals;
324 : : //! Unique set of ordered nodes at which symmetry BCs are set
325 : : std::set< std::size_t > m_symbcnodeset;
326 : : //! Streamable nodes at which symmetry BCs are set
327 : : std::vector< std::size_t > m_symbcnodes;
328 : : //! Streamable normals at nodes at which symmetry BCs are set
329 : : std::vector< tk::real > m_symbcnorms;
330 : : //! Unique set of ordered nodes at which farfield BCs are set
331 : : std::set< std::size_t > m_farbcnodeset;
332 : : //! Streamable nodes at which farfield BCs are set
333 : : std::vector< std::size_t > m_farbcnodes;
334 : : //! Streamable normals at nodes at which farfield BCs are set
335 : : std::vector< tk::real > m_farbcnorms;
336 : : //! Streamable surface integral nodes and normals * dA on surfaces
337 : : std::map< int, std::pair< std::vector< std::size_t >,
338 : : std::vector< tk::real > > > m_surfint;
339 : : //! Time step size for each mesh node
340 : : std::vector< tk::real > m_dtp;
341 : : //! Physical time for each mesh node
342 : : std::vector< tk::real > m_tp;
343 : : //! True in the last time step
344 : : int m_finished;
345 : : //! dt multiplier after flow no longer updated
346 : : tk::real m_freezeflow;
347 : : //! Freeze FCT limiter if 1, 0 FCT as usual
348 : : int m_fctfreeze;
349 : :
350 : : //! Access bound Discretization class pointer
351 : 303071 : Discretization* Disc() const {
352 : : Assert( m_disc[ thisIndex ].ckLocal() != nullptr, "ckLocal() null" );
353 : 606142 : return m_disc[ thisIndex ].ckLocal();
354 : : }
355 : :
356 : : //! Prepare boundary condition data structures
357 : : void setupBC();
358 : :
359 : : //! Compute local contributions to domain edge integrals
360 : : void domint();
361 : :
362 : : //! Compute chare-boundary edges
363 : : void bndEdges();
364 : :
365 : : //! Compute boundary point normals
366 : : void bndint();
367 : :
368 : : //! Combine own and communicated portions of the boundary point normals
369 : : void bnorm();
370 : :
371 : : //! Convert integrals into streamable data structures
372 : : void streamable();
373 : :
374 : : //! Generate superedge-groups for boundary-edge integrals
375 : : void bndsuped();
376 : :
377 : : //! Generate superedge-groups for domain-edge integral
378 : : void domsuped();
379 : :
380 : : //! Generate chare-boundary edge data structures for deactivation
381 : : void chbnded();
382 : :
383 : : //! Apply diffusion on active hull
384 : : void huldif();
385 : :
386 : : //! Output mesh and particle fields to files
387 : : void out();
388 : :
389 : : //! Output mesh-based fields to file
390 : : void writeFields( CkCallback cb );
391 : :
392 : : //! Combine own and communicated portions of the integrals
393 : : void merge();
394 : :
395 : : //! Compute righ-hand side vector of transport equations
396 : : void rhs();
397 : :
398 : : //! Continue with flux-corrected transport if enabled
399 : : void fct();
400 : :
401 : : //! Compute antidiffusive contributions: P+/-
402 : : void aec();
403 : :
404 : : //! Compute allowed limits, Q+/-
405 : : void alw();
406 : :
407 : : //! Compute limit coefficients
408 : : void lim();
409 : :
410 : : //! Advance systems of equations
411 : : void solve();
412 : :
413 : : //! Adjust node volumes along inactive neighbor chares
414 : : void deavol();
415 : :
416 : : //! Decide if edge is active
417 : : int active( std::size_t p,
418 : : std::size_t q,
419 : : const std::vector< uint64_t >& sys ) const;
420 : :
421 : : //! Decide whether to deactivate this chare
422 : : int dea( const std::vector< uint64_t >& sys ) const;
423 : :
424 : : //! Decide whether to teactivate a neighbor chare
425 : : void rea( const std::vector< uint64_t >& sys,
426 : : std::unordered_set< int >& req ) const;
427 : :
428 : : //! Deactivate regions
429 : : void deactivate();
430 : :
431 : : //! Compute deactivation status
432 : : void activate();
433 : :
434 : : //! Refine/derefine mesh
435 : : void refine();
436 : :
437 : : //! Compute time step size
438 : : void dt();
439 : :
440 : : //! Evaluate whether to save checkpoint/restart
441 : : void evalRestart();
442 : :
443 : : //! Apply boundary conditions
444 : : void BC( tk::Fields& u, tk::real t );
445 : :
446 : : //! Apply scalar source to solution
447 : : void src();
448 : : };
449 : :
450 : : } // inciter::
|
