Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/ALECG.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 ALECG for a PDE system with continuous Galerkin + ALE + RK
9 : : \details ALECG advances a system of partial differential equations (PDEs)
10 : : using a continuous Galerkin (CG) finite element (FE) spatial discretization
11 : : (using linear shapefunctions on tetrahedron elements) combined with a
12 : : Runge-Kutta (RK) time stepping scheme in the arbitrary Eulerian-Lagrangian
13 : : reference frame.
14 : :
15 : : There are a potentially large number of ALECG Charm++ chares created by
16 : : Transporter. Each ALECG gets a chunk of the full load (part of the mesh)
17 : : and does the same: initializes and advances a number of PDE systems in time.
18 : :
19 : : ALE time-stepping is performed in an unsplit fashion, as opposed to
20 : : Lagrange + remap. See also J. Waltz, N.R. Morgan, T.R. Canfield, M.R.J.
21 : : Charest, L.D. Risinger, J.G. Wohlbier, A three-dimensional finite element
22 : : arbitrary Lagrangian–Eulerian method for shock hydrodynamics on unstructured
23 : : grids, Computers & Fluids, 92: 172-187, 2014.
24 : :
25 : : The implementation uses the Charm++ runtime system and is fully
26 : : asynchronous, overlapping computation and communication. The algorithm
27 : : utilizes the structured dagger (SDAG) Charm++ functionality.
28 : : */
29 : : // *****************************************************************************
30 : : #ifndef ALECG_h
31 : : #define ALECG_h
32 : :
33 : : #include <vector>
34 : : #include <map>
35 : :
36 : : #include "Types.hpp"
37 : : #include "Fields.hpp"
38 : : #include "Table.hpp"
39 : : #include "DerivedData.hpp"
40 : : #include "NodeDiagnostics.hpp"
41 : : #include "Ghosts.hpp"
42 : :
43 : : #include "NoWarning/alecg.decl.h"
44 : :
45 : : namespace inciter {
46 : :
47 : : //! ALECG Charm++ chare array used to advance PDEs in time with ALECG+RK
48 : : class ALECG : public CBase_ALECG {
49 : :
50 : : public:
51 : : #if defined(__clang__)
52 : : #pragma clang diagnostic push
53 : : #pragma clang diagnostic ignored "-Wunused-parameter"
54 : : #pragma clang diagnostic ignored "-Wdeprecated-declarations"
55 : : #elif defined(STRICT_GNUC)
56 : : #pragma GCC diagnostic push
57 : : #pragma GCC diagnostic ignored "-Wunused-parameter"
58 : : #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
59 : : #elif defined(__INTEL_COMPILER)
60 : : #pragma warning( push )
61 : : #pragma warning( disable: 1478 )
62 : : #endif
63 : : // Include Charm++ SDAG code. See http://charm.cs.illinois.edu/manuals/html/
64 : : // charm++/manual.html, Sec. "Structured Control Flow: Structured Dagger".
65 : : ALECG_SDAG_CODE
66 : : #if defined(__clang__)
67 : : #pragma clang diagnostic pop
68 : : #elif defined(STRICT_GNUC)
69 : : #pragma GCC diagnostic pop
70 : : #elif defined(__INTEL_COMPILER)
71 : : #pragma warning( pop )
72 : : #endif
73 : :
74 : : //! Constructor
75 : : explicit ALECG( const CProxy_Discretization& disc,
76 : : const CProxy_Ghosts& ghostsproxy,
77 : : const std::map< int, std::vector< std::size_t > >& bface,
78 : : const std::map< int, std::vector< std::size_t > >& bnode,
79 : : const std::vector< std::size_t >& triinpoel );
80 : :
81 : : #if defined(__clang__)
82 : : #pragma clang diagnostic push
83 : : #pragma clang diagnostic ignored "-Wundefined-func-template"
84 : : #endif
85 : : //! Migrate constructor
86 : : // cppcheck-suppress uninitMemberVar
87 : 1624 : explicit ALECG( CkMigrateMessage* msg ) : CBase_ALECG( msg ) {}
88 : : #if defined(__clang__)
89 : : #pragma clang diagnostic pop
90 : : #endif
91 : :
92 : : //! Configure Charm++ custom reduction types initiated from this chare array
93 : : static void registerReducers();
94 : :
95 : : //! Return from migration
96 : : void ResumeFromSync() override;
97 : :
98 : : //! Size communication buffers (no-op)
99 : : void resizeComm() {}
100 : :
101 : : //! Setup node-neighborhood (no-op)
102 : : void nodeNeighSetup() {}
103 : :
104 : : //! (no-op)
105 : : void transferSol() {}
106 : :
107 : : //! Start setup for solution
108 : : void setup();
109 : :
110 : : //! Receive total box IC volume and set conditions in box
111 : : void box( tk::real v, const std::vector< tk::real >& blkvols );
112 : :
113 : : // Start time stepping
114 : : void start();
115 : :
116 : : //! Advance equations to next time step
117 : : void advance( tk::real newdt, tk::real );
118 : :
119 : : //! Compute left-hand side of transport equations
120 : : void lhs();
121 : :
122 : : //! Receive contributions to duual-face normals on chare boundaries
123 : : void comdfnorm(
124 : : const std::unordered_map< tk::UnsMesh::Edge,
125 : : std::array< tk::real, 3 >,
126 : : tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> >& dfnorm );
127 : :
128 : : //! Receive boundary point normals on chare-boundaries
129 : : void comnorm( const std::unordered_map< int,
130 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > >& innorm );
131 : :
132 : : //! Receive mesh block information for nodes on chare-boundaries
133 : : void comblk( const std::vector< std::size_t >& gid,
134 : : const std::vector< std::set< std::size_t > >& mb );
135 : :
136 : : //! Receive contributions to gradients on chare-boundaries
137 : : void comChBndGrad( const std::vector< std::size_t >& gid,
138 : : const std::vector< std::vector< tk::real > >& G );
139 : :
140 : : //! Receive contributions to right-hand side vector on chare-boundaries
141 : : void comrhs( const std::vector< std::size_t >& gid,
142 : : const std::vector< std::vector< tk::real > >& R );
143 : :
144 : : //! Update solution at the end of time step
145 : : void update( const tk::Fields& a );
146 : :
147 : : //! Optionally refine/derefine mesh
148 : : void refine( const std::vector< tk::real >& l2res );
149 : :
150 : : //! Receive new mesh from Refiner
151 : : void resizePostAMR(
152 : : const std::vector< std::size_t >& ginpoel,
153 : : const tk::UnsMesh::Chunk& chunk,
154 : : const tk::UnsMesh::Coords& coord,
155 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
156 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
157 : : const std::set< std::size_t >& removedNodes,
158 : : const std::unordered_map< std::size_t, std::size_t >& amrNodeMap,
159 : : const tk::NodeCommMap& nodeCommMap,
160 : : const std::map< int, std::vector< std::size_t > >& bface,
161 : : const std::map< int, std::vector< std::size_t > >& bnode,
162 : : const std::vector< std::size_t >& triinpoel,
163 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
164 : : elemblockid );
165 : :
166 : : //! Extract field output to file
167 : : void extractFieldOutput(
168 : : const std::vector< std::size_t >& /* ginpoel */,
169 : : const tk::UnsMesh::Chunk& /*chunk*/,
170 : : const tk::UnsMesh::Coords& /*coord*/,
171 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /* addedNodes */,
172 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
173 : : const tk::NodeCommMap& /*nodeCommMap*/,
174 : : const std::map< int, std::vector< std::size_t > >& /*bface*/,
175 : : const std::map< int, std::vector< std::size_t > >& /* bnode */,
176 : : const std::vector< std::size_t >& /*triinpoel*/,
177 : : CkCallback /*c*/ ) {}
178 : :
179 : : //! Const-ref access to current solution
180 : : //! \return Const-ref to current solution
181 : : const tk::Fields& solution() const { return m_u; }
182 : :
183 : : //! Start computing the mesh mesh velocity for ALE
184 : : void meshvelstart();
185 : :
186 : : //! Done with computing the mesh velocity for ALE mesh motion
187 : : void meshveldone();
188 : :
189 : : //! Resizing data sutrctures after mesh refinement has been completed
190 : : void resized();
191 : :
192 : : //! Evaluate whether to continue with next time step
193 : : void step();
194 : :
195 : : // Evaluate whether to do load balancing
196 : : void evalLB( int nrestart );
197 : :
198 : : //! Evaluate whether to continue with next time step stage
199 : : void stage();
200 : :
201 : : //! Continue to next time step
202 : : void next();
203 : :
204 : : /** @name Charm++ pack/unpack serializer member functions */
205 : : ///@{
206 : : //! \brief Pack/Unpack serialize member function
207 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
208 : 5112 : void pup( PUP::er &p ) override {
209 : 5112 : p | m_disc;
210 : 5112 : p | m_nsol;
211 : 5112 : p | m_ngrad;
212 : 5112 : p | m_nrhs;
213 : 5112 : p | m_nbnorm;
214 : 5112 : p | m_ndfnorm;
215 : 5112 : p | m_nmblk;
216 : 5112 : p | m_bnode;
217 : 5112 : p | m_bface;
218 : 5112 : p | m_triinpoel;
219 : 5112 : p | m_bndel;
220 : 5112 : p | m_dfnorm;
221 : 5112 : p | m_dfnormc;
222 : 5112 : p | m_dfn;
223 : 5112 : p | m_esup;
224 : 5112 : p | m_psup;
225 : 5112 : p | m_u;
226 : 5112 : p | m_un;
227 : 5112 : p | m_rhs;
228 : 5112 : p | m_rhsc;
229 : 5112 : p | m_chBndGrad;
230 : 5112 : p | m_dirbc;
231 : 5112 : p | m_chBndGradc;
232 : : p | m_diag;
233 : 5112 : p | m_bnorm;
234 : 5112 : p | m_bnormc;
235 : 5112 : p | m_symbcnodes;
236 : 5112 : p | m_farfieldbcnodes;
237 : 5112 : p | m_symbctri;
238 : 5112 : p | m_timedepbcnodes;
239 : 5112 : p | m_timedepbcFn;
240 : 5112 : p | m_stage;
241 : 5112 : p | m_boxnodes;
242 : 5112 : p | m_edgenode;
243 : 5112 : p | m_edgeid;
244 : 5112 : p | m_dtp;
245 : 5112 : p | m_tp;
246 : 5112 : p | m_finished;
247 : 5112 : p | m_newmesh;
248 : 5112 : p | m_refinedmesh;
249 : 5112 : p | m_nusermeshblk;
250 : 5112 : p | m_nodeblockid;
251 : 5112 : p | m_nodeblockidc;
252 : 5112 : }
253 : : //! \brief Pack/Unpack serialize operator|
254 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
255 : : //! \param[in,out] i ALECG object reference
256 : : friend void operator|( PUP::er& p, ALECG& i ) { i.pup(p); }
257 : : //@}
258 : :
259 : : private:
260 : : using ncomp_t = tk::ncomp_t;
261 : :
262 : : //! Discretization proxy
263 : : CProxy_Discretization m_disc;
264 : : //! Counter for high order solution vector nodes updated
265 : : std::size_t m_nsol;
266 : : //! Counter for nodal gradients updated
267 : : std::size_t m_ngrad;
268 : : //! Counter for right-hand side vector nodes updated
269 : : std::size_t m_nrhs;
270 : : //! Counter for receiving boundary point normals
271 : : std::size_t m_nbnorm;
272 : : //! Counter for receiving dual-face normals on chare-boundary edges
273 : : std::size_t m_ndfnorm;
274 : : //! Counter for mesh block nodes updated
275 : : std::size_t m_nmblk;
276 : : //! Boundary node lists mapped to side set ids used in the input file
277 : : std::map< int, std::vector< std::size_t > > m_bnode;
278 : : //! Boundary face lists mapped to side set ids used in the input file
279 : : std::map< int, std::vector< std::size_t > > m_bface;
280 : : //! Boundary triangle face connecitivity where BCs are set by user
281 : : std::vector< std::size_t > m_triinpoel;
282 : : //! Elements along mesh boundary
283 : : std::vector< std::size_t > m_bndel;
284 : : //! Dual-face normals along edges
285 : : std::unordered_map< tk::UnsMesh::Edge, std::array< tk::real, 3 >,
286 : : tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> > m_dfnorm;
287 : : //! Receive buffer for dual-face normals along chare-boundary edges
288 : : std::unordered_map< tk::UnsMesh::Edge, std::array< tk::real, 3 >,
289 : : tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> > m_dfnormc;
290 : : //! Streamable dual-face normals
291 : : std::vector< tk::real > m_dfn;
292 : : //! El;ements surrounding points
293 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > > m_esup;
294 : : //! Points surrounding points
295 : : std::pair< std::vector< std::size_t >, std::vector< std::size_t > > m_psup;
296 : : //! Unknown/solution vector at mesh nodes
297 : : tk::Fields m_u;
298 : : //! Unknown/solution vector at mesh nodes at previous time
299 : : tk::Fields m_un;
300 : : //! Right-hand side vector (for the high order system)
301 : : tk::Fields m_rhs;
302 : : //! Receive buffer for communication of the right hand side
303 : : //! \details Key: global node id, value: rhs for all scalar components per
304 : : //! node.
305 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_rhsc;
306 : : //! Nodal gradients at chare-boundary nodes
307 : : tk::Fields m_chBndGrad;
308 : : //! Boundary conditions evaluated and assigned to local mesh node IDs
309 : : //! \details Vector of pairs of bool and boundary condition value associated
310 : : //! to local mesh node IDs at which the user has set Dirichlet boundary
311 : : //! conditions for all PDEs integrated. The bool indicates whether the BC
312 : : //! is set at the node for that component the if true, the real value is
313 : : //! the increment (from t to dt) in the BC specified for a component.
314 : : std::unordered_map< std::size_t,
315 : : std::vector< std::pair< bool, tk::real > > > m_dirbc;
316 : : //! Receive buffer for communication of the nodal gradients
317 : : //! \details Key: chare id, value: gradients for all scalar components per
318 : : //! node
319 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_chBndGradc;
320 : : //! Diagnostics object
321 : : NodeDiagnostics m_diag;
322 : : //! Face normals in boundary points associated to side sets
323 : : //! \details Key: local node id, value: unit normal and inverse distance
324 : : //! square between face centroids and points, outer key: side set id
325 : : std::unordered_map< int,
326 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > > m_bnorm;
327 : : //! \brief Receive buffer for communication of the boundary point normals
328 : : //! associated to side sets
329 : : //! \details Key: global node id, value: normals (first 3 components),
330 : : //! inverse distance squared (4th component), outer key, side set id
331 : : std::unordered_map< int,
332 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > > m_bnormc;
333 : : //! Unique set of nodes at which symmetry BCs are set
334 : : std::unordered_set< std::size_t > m_symbcnodes;
335 : : //! Unique set of nodes at which farfield BCs are set
336 : : std::unordered_set< std::size_t > m_farfieldbcnodes;
337 : : //! Vector with 1 at symmetry BC boundary triangles
338 : : std::vector< int > m_symbctri;
339 : : //! \brief Unique set of nodes at which time dependent BCs are set
340 : : // for each time dependent BC
341 : : std::vector< std::unordered_set< std::size_t > > m_timedepbcnodes;
342 : : //! \brief User defined discrete function of time used in the time dependent
343 : : // BCs associated with (index in vector) the number of distinct time
344 : : // dependent BCs specified. This index is the same as the index in
345 : : // m_timedepbcnodes.
346 : : std::vector< tk::Table<5> > m_timedepbcFn;
347 : : //! Runge-Kutta stage counter
348 : : std::size_t m_stage;
349 : : //! Mesh node ids at which user-defined box ICs are defined (multiple boxes)
350 : : std::vector< std::unordered_set< std::size_t > > m_boxnodes;
351 : : //! Local node IDs of edges
352 : : std::vector< std::size_t > m_edgenode;
353 : : //! Edge ids in the order of access
354 : : std::vector< std::size_t > m_edgeid;
355 : : //! Time step size for each mesh node
356 : : std::vector< tk::real > m_dtp;
357 : : //! Physical time for each mesh node
358 : : std::vector< tk::real > m_tp;
359 : : //! True in the last time step
360 : : int m_finished;
361 : : //! \brief State indicating the reason we are recomputing the normals.
362 : : // 0: after ALE; 1: after AMR
363 : : int m_newmesh;
364 : : //! State indicating if the mesh has been refined by dtref
365 : : int m_refinedmesh;
366 : : //! Number of mesh-blocks with user-defined ICs
367 : : std::size_t m_nusermeshblk;
368 : : //! Local node ids associated with mesh block ids
369 : : std::unordered_map< std::size_t, std::set< std::size_t > > m_nodeblockid;
370 : : //! Receive buffer for communication of the mesh block ids
371 : : //! \details Key: mesh block id, value: set of global node ids for nodes
372 : : //! in this mesh block.
373 : : std::unordered_map< std::size_t, std::set< std::size_t > > m_nodeblockidc;
374 : :
375 : : //! Access bound Discretization class pointer
376 : 9257803 : Discretization* Disc() const {
377 : : Assert( m_disc[ thisIndex ].ckLocal() != nullptr, "ckLocal() null" );
378 : 18515606 : return m_disc[ thisIndex ].ckLocal();
379 : : }
380 : :
381 : : //! Compute normal of dual-mesh associated to edge
382 : : std::array< tk::real, 3 >
383 : : edfnorm( const tk::UnsMesh::Edge& edge,
384 : : const std::unordered_map< tk::UnsMesh::Edge,
385 : : std::vector< std::size_t >,
386 : : tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> >& esued ) const;
387 : :
388 : : //! Compute chare-boundary edges
389 : : void bndEdges();
390 : :
391 : : //! Start (re-)computing boundare point-, and dual-face normals
392 : : void norm();
393 : :
394 : : //! Compute dual-face normals associated to edges
395 : : void dfnorm();
396 : :
397 : : //! Compute boundary point normals
398 : : void
399 : : bnorm( const std::unordered_map< int,
400 : : std::unordered_set< std::size_t > >& bcnodes );
401 : :
402 : : //! \brief Finish computing dual-face and boundary point normals and apply
403 : : //! boundary conditions on the initial conditions
404 : : void normfinal();
405 : :
406 : : //! Continue setup for solution, after communication for mesh blocks
407 : : void continueSetup();
408 : :
409 : : //! Output mesh and particle fields to files
410 : : void out();
411 : :
412 : : //! Output mesh-based fields to file
413 : : void writeFields( CkCallback c );
414 : :
415 : : //! Combine own and communicated contributions to normals
416 : : void mergelhs();
417 : :
418 : : //! Compute gradients
419 : : void chBndGrad();
420 : :
421 : : //! Compute righ-hand side vector of transport equations
422 : : void rhs();
423 : :
424 : : //! Advance systems of equations
425 : : void solve();
426 : :
427 : : //! Compute time step size
428 : : void dt();
429 : :
430 : : //! Transfer solution to other solver and mesh if coupled
431 : : void transfer();
432 : :
433 : : //! Evaluate whether to save checkpoint/restart
434 : : void evalRestart();
435 : :
436 : : //! Query/update boundary-conditions-related data structures from user input
437 : : void queryBnd();
438 : :
439 : : //! Apply boundary conditions
440 : : void BC();
441 : :
442 : : //! Multiply solution with mesh volume
443 : : void volumetric( tk::Fields& u, const std::vector< tk::real >& v );
444 : :
445 : : //! Divide solution with mesh volume
446 : : void conserved( tk::Fields& u, const std::vector< tk::real >& v );
447 : :
448 : : //! Continue after computing the new mesh velocity for ALE
449 : : void ale();
450 : : };
451 : :
452 : : } // inciter::
453 : :
454 : : #endif // ALECG_h
|
