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