Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/DiagCG.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 DiagCG for a PDE system with continuous Galerkin without a matrix
9 : : \details DiagCG advances a system of partial differential equations (PDEs)
10 : : using continuous Galerkin (CG) finite element (FE) spatial discretization
11 : : (using linear shapefunctions on tetrahedron elements) combined with a time
12 : : stepping scheme that is equivalent to the Lax-Wendroff (LW) scheme within
13 : : the unstructured-mesh FE context and treats discontinuities with
14 : : flux-corrected transport (FCT). The left-hand side (lumped-mass) matrix
15 : : is diagonal thus this scheme does not use a matrix-based linear solver.
16 : :
17 : : There are a potentially large number of DiagCG Charm++ chares created by
18 : : Transporter. Each DiagCG gets a chunk of the full load (part of the mesh)
19 : : and does the same: initializes and advances a number of PDE systems in time.
20 : :
21 : : The implementation uses the Charm++ runtime system and is fully
22 : : asynchronous, overlapping computation and communication. The algorithm
23 : : utilizes the structured dagger (SDAG) Charm++ functionality. The high-level
24 : : overview of the algorithm structure and how it interfaces with Charm++ is
25 : : discussed in the Charm++ interface file src/Inciter/diagcg.ci.
26 : : */
27 : : // *****************************************************************************
28 : : #ifndef DiagCG_h
29 : : #define DiagCG_h
30 : :
31 : : #include <vector>
32 : : #include <map>
33 : : #include <unordered_set>
34 : :
35 : : #include "Types.hpp"
36 : : #include "Fields.hpp"
37 : : #include "DerivedData.hpp"
38 : : #include "FluxCorrector.hpp"
39 : : #include "NodeDiagnostics.hpp"
40 : : #include "CommMap.hpp"
41 : : #include "Inciter/InputDeck/InputDeck.hpp"
42 : :
43 : : #include "NoWarning/diagcg.decl.h"
44 : :
45 : : namespace inciter {
46 : :
47 : : extern ctr::InputDeck g_inputdeck;
48 : :
49 : : //! DiagCG Charm++ chare array used to advance PDEs in time with DiagCG+LW+FCT
50 : : class DiagCG : public CBase_DiagCG {
51 : :
52 : : public:
53 : : #if defined(__clang__)
54 : : #pragma clang diagnostic push
55 : : #pragma clang diagnostic ignored "-Wunused-parameter"
56 : : #pragma clang diagnostic ignored "-Wdeprecated-declarations"
57 : : #elif defined(STRICT_GNUC)
58 : : #pragma GCC diagnostic push
59 : : #pragma GCC diagnostic ignored "-Wunused-parameter"
60 : : #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
61 : : #elif defined(__INTEL_COMPILER)
62 : : #pragma warning( push )
63 : : #pragma warning( disable: 1478 )
64 : : #endif
65 : : // Include Charm++ SDAG code. See http://charm.cs.illinois.edu/manuals/html/
66 : : // charm++/manual.html, Sec. "Structured Control Flow: Structured Dagger".
67 : : DiagCG_SDAG_CODE
68 : : #if defined(__clang__)
69 : : #pragma clang diagnostic pop
70 : : #elif defined(STRICT_GNUC)
71 : : #pragma GCC diagnostic pop
72 : : #elif defined(__INTEL_COMPILER)
73 : : #pragma warning( pop )
74 : : #endif
75 : :
76 : : //! Constructor
77 : : explicit DiagCG( const CProxy_Discretization& disc,
78 : : const std::map< int, std::vector< std::size_t > >& bface,
79 : : const std::map< int, std::vector< std::size_t > >& bnode,
80 : : const std::vector< std::size_t >& triinpoel );
81 : :
82 : : #if defined(__clang__)
83 : : #pragma clang diagnostic push
84 : : #pragma clang diagnostic ignored "-Wundefined-func-template"
85 : : #endif
86 : : //! Migrate constructor
87 : : // cppcheck-suppress uninitMemberVar
88 : 1896 : explicit DiagCG( CkMigrateMessage* msg ) : CBase_DiagCG( msg ) {}
89 : : #if defined(__clang__)
90 : : #pragma clang diagnostic pop
91 : : #endif
92 : :
93 : : //! Configure Charm++ custom reduction types initiated from this chare array
94 : : static void registerReducers();
95 : :
96 : : //! Return from migration
97 : : void ResumeFromSync() override;
98 : :
99 : : //! Size communication buffers (no-op)
100 : : void resizeComm() {}
101 : :
102 : : //! Setup node-neighborhood (no-op)
103 : : void nodeNeighSetup() {}
104 : :
105 : : //! Setup: query boundary conditions, output mesh, etc.
106 : : void setup();
107 : :
108 : : //! Receive total box IC volume and set conditions in box
109 : : void box( tk::real v );
110 : :
111 : : // Initially compute left hand side diagonal matrix
112 : : void init();
113 : :
114 : : //! Advance equations to next time step
115 : : void advance( tk::real newdt );
116 : :
117 : : //! Compute left-hand side of transport equations
118 : : void lhs();
119 : :
120 : : //! Receive boundary point normals on chare-boundaries
121 : : void comnorm( const std::unordered_map< int,
122 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > >& innorm );
123 : :
124 : : //! Receive contributions to left-hand side matrix on chare-boundaries
125 : : void comlhs( const std::vector< std::size_t >& gid,
126 : : const std::vector< std::vector< tk::real > >& L );
127 : :
128 : : //! Receive contributions to right-hand side vector on chare-boundaries
129 : : void comrhs( const std::vector< std::size_t >& gid,
130 : : const std::vector< std::vector< tk::real > >& R,
131 : : const std::vector< std::vector< tk::real > >& D );
132 : :
133 : : //! Update solution at the end of time step
134 : : void update( const tk::Fields& a, tk::Fields&& dul );
135 : :
136 : : //! Optionally refine/derefine mesh
137 : : void refine( const std::vector< tk::real >& l2res );
138 : :
139 : : //! Receive new mesh from Refiner
140 : : void resizePostAMR(
141 : : const std::vector< std::size_t >& ginpoel,
142 : : const tk::UnsMesh::Chunk& chunk,
143 : : const tk::UnsMesh::Coords& coord,
144 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
145 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
146 : : const std::set< std::size_t >& removedNodes,
147 : : const tk::NodeCommMap& nodeCommMap,
148 : : const std::map< int, std::vector< std::size_t > >& /* bface */,
149 : : const std::map< int, std::vector< std::size_t > >& bnode,
150 : : const std::vector< std::size_t >& /* triinpoel */ );
151 : :
152 : : //! Extract field output to file
153 : : void extractFieldOutput(
154 : : const std::vector< std::size_t >& /* ginpoel */,
155 : : const tk::UnsMesh::Chunk& /*chunk*/,
156 : : const tk::UnsMesh::Coords& /*coord*/,
157 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /* addedNodes */,
158 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
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 : : CkCallback /*c*/ ) {}
164 : :
165 : : //! Const-ref access to current solution
166 : : //! \return Const-ref to current solution
167 : : const tk::Fields& solution() const { return m_u; }
168 : :
169 : : //! Resizing data sutrctures after mesh refinement has been completed
170 : : void resized();
171 : :
172 : : //! Evaluate whether to continue with next time step
173 : : void step();
174 : :
175 : : // Evaluate whether to do load balancing
176 : : void evalLB( int nrestart );
177 : :
178 : : //! Continue to next time step
179 : : void next();
180 : :
181 : : /** @name Charm++ pack/unpack serializer member functions */
182 : : ///@{
183 : : //! \brief Pack/Unpack serialize member function
184 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
185 : 5688 : void pup( PUP::er &p ) override {
186 : : p | m_disc;
187 : 5688 : p | m_initial;
188 : 5688 : p | m_nlhs;
189 : 5688 : p | m_nrhs;
190 : 5688 : p | m_nnorm;
191 : 5688 : p | m_bnode;
192 : 5688 : p | m_bface;
193 : 5688 : p | m_triinpoel;
194 : 5688 : p | m_u;
195 : 5688 : p | m_ul;
196 : 5688 : p | m_du;
197 : 5688 : p | m_ue;
198 : 5688 : p | m_lhs;
199 : 5688 : p | m_rhs;
200 : 5688 : p | m_bcdir;
201 : 5688 : p | m_lhsc;
202 : 5688 : p | m_rhsc;
203 : 5688 : p | m_difc;
204 : 5688 : p | m_vol;
205 : 5688 : p | m_bnorm;
206 : 5688 : p | m_bnormc;
207 : 5688 : p | m_symbcnodemap;
208 : 5688 : p | m_symbcnodes;
209 : 5688 : p | m_farfieldbcnodes;
210 : : p | m_diag;
211 : 5688 : p | m_boxnodes;
212 : 5688 : p | m_dtp;
213 : 5688 : p | m_tp;
214 : 5688 : }
215 : : //! \brief Pack/Unpack serialize operator|
216 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
217 : : //! \param[in,out] i DiagCG object reference
218 : : friend void operator|( PUP::er& p, DiagCG& i ) { i.pup(p); }
219 : : //@}
220 : :
221 : : private:
222 : : using ncomp_t = kw::ncomp::info::expect::type;
223 : :
224 : : //! Discretization proxy
225 : : CProxy_Discretization m_disc;
226 : : //! 1 if starting time stepping, 0 if during time stepping
227 : : std::size_t m_initial;
228 : : //! Counter for left-hand side matrix (vector) nodes updated
229 : : std::size_t m_nlhs;
230 : : //! Counter for right-hand side vector nodes updated
231 : : std::size_t m_nrhs;
232 : : //! Counter for receiving boundary point normals
233 : : std::size_t m_nnorm;
234 : : //! Boundary node lists mapped to side set ids
235 : : std::map< int, std::vector< std::size_t > > m_bnode;
236 : : //! Boundary faces side-set information
237 : : std::map< int, std::vector< std::size_t > > m_bface;
238 : : //! Triangle face connecitivity
239 : : std::vector< std::size_t > m_triinpoel;
240 : : //! Unknown/solution vector at mesh nodes
241 : : tk::Fields m_u;
242 : : //! Unknown/solution vector at mesh nodes (low orderd)
243 : : tk::Fields m_ul;
244 : : //! Unknown/solution vector increment (high order)
245 : : tk::Fields m_du;
246 : : //! Unknown/solution vector at mesh cells
247 : : tk::Fields m_ue;
248 : : //! Lumped lhs mass matrix
249 : : tk::Fields m_lhs;
250 : : //! Right-hand side vector (for the high order system)
251 : : tk::Fields m_rhs;
252 : : //! Boundary conditions evaluated and assigned to local mesh node IDs
253 : : //! \details Vector of pairs of bool and boundary condition value associated
254 : : //! to local mesh node IDs at which the user has set Dirichlet boundary
255 : : //! conditions for all PDEs integrated. The bool indicates whether the BC
256 : : //! is set at the node for that component the if true, the real value is
257 : : //! the increment (from t to dt) in the BC specified for a component.
258 : : std::unordered_map< std::size_t,
259 : : std::vector< std::pair< bool, tk::real > > > m_bcdir;
260 : : //! Receive buffer for communication of the left hand side
261 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_lhsc;
262 : : //! Receive buffer for communication of the right hand side
263 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_rhsc;
264 : : //! Receive buffer for communication of mass diffusion on the hand side
265 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_difc;
266 : : //! Total mesh volume
267 : : tk::real m_vol;
268 : : //! Face normals in boundary points associated to side sets
269 : : //! \details Key: local node id, value: unit normal and inverse distance
270 : : //! square between face centroids and points, outer key: side set id
271 : : std::unordered_map< int,
272 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > > m_bnorm;
273 : : //! \brief Receive buffer for communication of the boundary point normals
274 : : //! associated to side sets
275 : : //! \details Key: global node id, value: normals (first 3 components),
276 : : //! inverse distance squared (4th component), outer key, side set id
277 : : std::unordered_map< int,
278 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > > m_bnormc;
279 : : //! Unique set of nodes at which symmetry BCs are set for side sets
280 : : std::unordered_map< int, std::unordered_set< std::size_t > > m_symbcnodemap;
281 : : //! Unique set of nodes at which symmetry BCs are set
282 : : std::unordered_set< std::size_t > m_symbcnodes;
283 : : //! Unique set of nodes at which farfield BCs are set
284 : : std::unordered_set< std::size_t > m_farfieldbcnodes;
285 : : //! Diagnostics object
286 : : NodeDiagnostics m_diag;
287 : : //! Mesh node ids at which user-defined box ICs are defined (multiple boxes)
288 : : std::vector< std::unordered_set< std::size_t > > m_boxnodes;
289 : : //! Time step size for each mesh node
290 : : std::vector< tk::real > m_dtp;
291 : : //! Physical time for each mesh node
292 : : std::vector< tk::real > m_tp;
293 : : //! True in the last time step
294 : : int m_finished;
295 : :
296 : : //! Access bound Discretization class pointer
297 : 469542 : Discretization* Disc() const {
298 : : Assert( m_disc[ thisIndex ].ckLocal() != nullptr, "ckLocal() null" );
299 : 939084 : return m_disc[ thisIndex ].ckLocal();
300 : : }
301 : :
302 : : //! Compute boundary point normals
303 : : void bnorm( const std::unordered_map< int,
304 : : std::unordered_set< std::size_t > >& bcnodes );
305 : :
306 : : //! Finish setting up communication maps (norms, etc.)
307 : : void normfinal();
308 : :
309 : : //! Output mesh fields to files
310 : : void out();
311 : :
312 : : //! Output mesh-based fields to file
313 : : void writeFields( CkCallback c ) const;
314 : :
315 : : //! The own and communication portion of the left-hand side is complete
316 : : void lhsdone();
317 : :
318 : : //! Combine own and communicated contributions to left hand side
319 : : void lhsmerge();
320 : :
321 : : //! Compute righ-hand side vector of transport equations
322 : : void rhs();
323 : :
324 : : //! Start time stepping
325 : : void start();
326 : :
327 : : //! Solve low and high order diagonal systems
328 : : void solve( tk::Fields& dif );
329 : :
330 : : //! Compute time step size
331 : : void dt();
332 : :
333 : : //! Evaluate whether to save checkpoint/restart
334 : : void evalRestart();
335 : : };
336 : :
337 : : } // inciter::
338 : :
339 : : #endif // DiagCG_h
|
