Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/LohCG.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 LohCG: Artificial compressibility solver for incompressible flow
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 "PUPUtil.hpp"
24 : :
25 : : #include "NoWarning/lohcg.decl.h"
26 : :
27 : : namespace inciter {
28 : :
29 : : //! LohCG Charm++ chare array used to advance PDEs in time with LohCG
30 : : class LohCG : public CBase_LohCG {
31 : :
32 : : public:
33 : : #if defined(__clang__)
34 : : #pragma clang diagnostic push
35 : : #pragma clang diagnostic ignored "-Wunused-parameter"
36 : : #pragma clang diagnostic ignored "-Wdeprecated-declarations"
37 : : #elif defined(STRICT_GNUC)
38 : : #pragma GCC diagnostic push
39 : : #pragma GCC diagnostic ignored "-Wunused-parameter"
40 : : #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
41 : : #elif defined(__INTEL_COMPILER)
42 : : #pragma warning( push )
43 : : #pragma warning( disable: 1478 )
44 : : #endif
45 : : // Include Charm++ SDAG code. See http://charm.cs.illinois.edu/manuals/html/
46 : : // charm++/manual.html, Sec. "Structured Control Flow: Structured Dagger".
47 : : LohCG_SDAG_CODE
48 : : #if defined(__clang__)
49 : : #pragma clang diagnostic pop
50 : : #elif defined(STRICT_GNUC)
51 : : #pragma GCC diagnostic pop
52 : : #elif defined(__INTEL_COMPILER)
53 : : #pragma warning( pop )
54 : : #endif
55 : :
56 : : //! Constructor
57 : : explicit LohCG( const CProxy_Discretization& disc,
58 : : const tk::CProxy_ConjugateGradients& cgpre,
59 : : const std::map< int, std::vector< std::size_t > >& bface,
60 : : const std::map< int, std::vector< std::size_t > >& bnode,
61 : : const std::vector< std::size_t >& triinpoel );
62 : :
63 : : #if defined(__clang__)
64 : : #pragma clang diagnostic push
65 : : #pragma clang diagnostic ignored "-Wundefined-func-template"
66 : : #endif
67 : : //! Migrate constructor
68 : : // cppcheck-suppress uninitMemberVar
69 : 653 : explicit LohCG( CkMigrateMessage* m ) : CBase_LohCG( m ) {}
70 : : #if defined(__clang__)
71 : : #pragma clang diagnostic pop
72 : : #endif
73 : :
74 : : //! Configure Charm++ custom reduction types initiated from this chare array
75 : : static void registerReducers();
76 : :
77 : : //! Return from migration
78 : : void ResumeFromSync() override;
79 : :
80 : : //! Start setup for solution
81 : : void setup( tk::real v );
82 : :
83 : : //! Initialize Poisson solve
84 : : void pinit();
85 : :
86 : : //! Solve Poisson equation
87 : : void psolve();
88 : :
89 : : //! Continue after Poisson solve
90 : : void psolved();
91 : :
92 : : // Start time stepping
93 : : void start();
94 : :
95 : : //! Advance equations to next time step
96 : : void advance( tk::real newdt );
97 : :
98 : : //! Evaluate diagnostics
99 : : void diag();
100 : :
101 : : //! Start (re-)computing domain and boundary integrals
102 : : void feop();
103 : :
104 : : //! Receive contributions to boundary point normals on chare-boundaries
105 : : void comnorm( const std::unordered_map< int,
106 : : std::unordered_map< std::size_t, std::array<tk::real,4> > >& inbnd );
107 : :
108 : : //! Receive contributions to velocity gradients
109 : : void comvgrad( const std::unordered_map< std::size_t,
110 : : std::vector< tk::real > >& ingrad );
111 : :
112 : : //! Receive contributions to momentum flux on chare-boundaries
113 : : void comflux( const std::unordered_map< std::size_t,
114 : : std::vector< tk::real > >& influx );
115 : :
116 : : //! Receive contributions to conjugate gradients solution gradient
117 : : void comsgrad( const std::unordered_map< std::size_t,
118 : : std::vector< tk::real > >& ingrad );
119 : :
120 : : //! Receive contributions to gradient
121 : : void comgrad( const std::unordered_map< std::size_t,
122 : : std::vector< tk::real > >& ingrad );
123 : :
124 : : //! Receive contributions to right-hand side vector on chare-boundaries
125 : : void comrhs( const std::unordered_map< std::size_t,
126 : : std::vector< tk::real > >& inrhs );
127 : :
128 : : //! Receive contributions to velocity divergence on chare-boundaries
129 : : void comdiv( const std::unordered_map< std::size_t, tk::real >& indiv );
130 : :
131 : : //! Evaluate residuals
132 : : void evalres( const std::vector< tk::real >& l2res );
133 : :
134 : : //! Receive new mesh from Refiner
135 : : void resizePostAMR(
136 : : const std::vector< std::size_t >& ginpoel,
137 : : const tk::UnsMesh::Chunk& chunk,
138 : : const tk::UnsMesh::Coords& coord,
139 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
140 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
141 : : const std::set< std::size_t >& removedNodes,
142 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
143 : : nodeCommMap,
144 : : const std::map< int, std::vector< std::size_t > >& bface,
145 : : const std::map< int, std::vector< std::size_t > >& bnode,
146 : : const std::vector< std::size_t >& triinpoel );
147 : :
148 : : //! Const-ref access to current solution
149 : : //! \return Const-ref to current solution
150 : 0 : const tk::Fields& solution() const { return m_u; }
151 : :
152 : : //! Compute integral quantities for output
153 : : void integrals();
154 : :
155 : : //! Compute recent conjugate gradients solution gradient
156 : : void sgrad();
157 : :
158 : : //! Evaluate whether to continue with next time step
159 : : void step();
160 : :
161 : : // Evaluate whether to do load balancing
162 : : void evalLB( int nrestart );
163 : :
164 : : /** @name Charm++ pack/unpack serializer member functions */
165 : : ///@{
166 : : //! \brief Pack/Unpack serialize member function
167 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
168 : 1987 : void pup( PUP::er &p ) override {
169 : 1987 : p | m_disc;
170 : 1987 : p | m_cgpre;
171 : 1987 : p | m_nrhs;
172 : 1987 : p | m_nnorm;
173 : 1987 : p | m_ngrad;
174 : 1987 : p | m_nsgrad;
175 : 1987 : p | m_nvgrad;
176 : 1987 : p | m_nflux;
177 : 1987 : p | m_ndiv;
178 : 1987 : p | m_nbpint;
179 : 1987 : p | m_np;
180 : 1987 : p | m_bnode;
181 : 1987 : p | m_bface;
182 : 1987 : p | m_triinpoel;
183 : 1987 : p | m_u;
184 : 1987 : p | m_un;
185 : 1987 : p | m_grad;
186 : 1987 : p | m_gradc;
187 : 1987 : p | m_vgrad;
188 : 1987 : p | m_vgradc;
189 : 1987 : p | m_flux;
190 : 1987 : p | m_fluxc;
191 : 1987 : p | m_div;
192 : 1987 : p | m_divc;
193 : 1987 : p | m_sgrad;
194 : 1987 : p | m_sgradc;
195 : 1987 : p | m_rhs;
196 : 1987 : p | m_rhsc;
197 : 1987 : p | m_diag;
198 : 1987 : p | m_bnorm;
199 : 1987 : p | m_bnormc;
200 : 1987 : p | m_bndpoinint;
201 : 1987 : p | m_domedgeint;
202 : 1987 : p | m_bpint;
203 : 1987 : p | m_dsupedge;
204 : 1987 : p | m_dsupint;
205 : 1987 : p | m_dirbcmask;
206 : 1987 : p | m_dirbcval;
207 : 1987 : p | m_dirbcmaskp;
208 : 1987 : p | m_dirbcvalp;
209 : 1987 : p | m_symbcnodes;
210 : 1987 : p | m_symbcnorms;
211 : 1987 : p | m_noslipbcnodes;
212 : 1987 : p | m_surfint;
213 : 1987 : p | m_stage;
214 : 1987 : p | m_finished;
215 : 1987 : p | m_rkcoef;
216 : 1987 : p | m_timer;
217 : 1987 : }
218 : : //! \brief Pack/Unpack serialize operator|
219 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
220 : : //! \param[in,out] i LohCG object reference
221 : : friend void operator|( PUP::er& p, LohCG& i ) { i.pup(p); }
222 : : //@}
223 : :
224 : : private:
225 : : //! Discretization proxy
226 : : CProxy_Discretization m_disc;
227 : : //! Conjugate Gradients Charm++ proxy for pressure solve
228 : : tk::CProxy_ConjugateGradients m_cgpre;
229 : : //! Counter for right-hand side vector nodes updated
230 : : std::size_t m_nrhs;
231 : : //! Counter for receiving boundary point normals
232 : : std::size_t m_nnorm;
233 : : //! Counter for receiving gradient
234 : : std::size_t m_ngrad;
235 : : //! Counter for receiving conjugrate gradient solution gradient
236 : : std::size_t m_nsgrad;
237 : : //! Counter for receiving velocity gradient
238 : : std::size_t m_nvgrad;
239 : : //! Counter for receiving momentum flux
240 : : std::size_t m_nflux;
241 : : //! Counter for receiving boundary velocity divergences
242 : : std::size_t m_ndiv;
243 : : //! Counter for receiving boundary point integrals
244 : : std::size_t m_nbpint;
245 : : //! Count number of Poisson solves during setup
246 : : std::size_t m_np;
247 : : //! Boundary node lists mapped to side set ids used in the input file
248 : : std::map< int, std::vector< std::size_t > > m_bnode;
249 : : //! Boundary face lists mapped to side set ids used in the input file
250 : : std::map< int, std::vector< std::size_t > > m_bface;
251 : : //! Boundary triangle face connecitivity where BCs are set by user
252 : : std::vector< std::size_t > m_triinpoel;
253 : : //! Unknown/solution vector at mesh nodes
254 : : tk::Fields m_u;
255 : : //! Unknown/solution vector at mesh nodes at previous time step
256 : : tk::Fields m_un;
257 : : //! Gradient in mesh nodes
258 : : tk::Fields m_grad;
259 : : //! Gradient receive buffer
260 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_gradc;
261 : : //! Velocity gradient in mesh nodes
262 : : tk::Fields m_vgrad;
263 : : //! Velocity gradient receive buffer
264 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_vgradc;
265 : : //! Momentum flux in mesh nodes
266 : : tk::Fields m_flux;
267 : : //! Momentum flux receive buffer
268 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_fluxc;
269 : : //! Velocity divergence
270 : : std::vector< tk::real > m_div;
271 : : //! Receive buffer for communication of the velocity divergence
272 : : //! \details Key: global node id, value: velocity divergence
273 : : std::unordered_map< std::size_t, tk::real > m_divc;
274 : : //! Conjugate gradient solution gradient in mesh nodes
275 : : tk::Fields m_sgrad;
276 : : //! Conjugate gradient solution gradient receive buffer
277 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_sgradc;
278 : : //! Right-hand side vector (for the high order system)
279 : : tk::Fields m_rhs;
280 : : //! Receive buffer for communication of the right hand side
281 : : //! \details Key: global node id, value: rhs for all scalar components per
282 : : //! node.
283 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_rhsc;
284 : : //! Diagnostics object
285 : : NodeDiagnostics m_diag;
286 : : //! Boundary point normals
287 : : //! \details Outer key: side set id. Inner key: global node id of boundary
288 : : //! point, value: weighted normal vector, inverse distance square.
289 : : std::unordered_map< int,
290 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > > m_bnorm;
291 : : //! Boundary point normals receive buffer
292 : : //! \details Outer key: side set id. Inner key: global node id of boundary
293 : : //! point, value: weighted normals and inverse distance square.
294 : : decltype(m_bnorm) m_bnormc;
295 : : //! Boundary point integrals
296 : : //! \details Key: global node id of boundary point, value: boundary point
297 : : //! integral contributions.
298 : : std::unordered_map< std::size_t, std::array<tk::real,3> > m_bndpoinint;
299 : : //! Domain edge integrals
300 : : std::unordered_map< tk::UnsMesh::Edge, std::array< tk::real, 4 >,
301 : : tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> > m_domedgeint;
302 : : //! Streamable boundary point integrals
303 : : std::vector< tk::real > m_bpint;
304 : : //! Superedge (tet, face, edge) end points with local ids for domain edges
305 : : std::array< std::vector< std::size_t >, 3 > m_dsupedge;
306 : : //! Superedge (tet, face, edge) domain edge integrals
307 : : std::array< std::vector< tk::real >, 3 > m_dsupint;
308 : : //! Nodes and their Dirichlet BC masks
309 : : std::vector< std::size_t > m_dirbcmask;
310 : : //! Nodes and their Dirichlet BC values
311 : : std::vector< double > m_dirbcval;
312 : : //! Nodes and their pressure Dirichlet BC masks
313 : : std::vector< std::size_t > m_dirbcmaskp;
314 : : //! Nodes and their pressure Dirichlet BC values
315 : : std::vector< double > m_dirbcvalp;
316 : : //! Streamable nodes at which symmetry BCs are set
317 : : std::vector< std::size_t > m_symbcnodes;
318 : : //! Streamable normals at nodes at which symmetry BCs are set
319 : : std::vector< tk::real > m_symbcnorms;
320 : : //! Streamable nodes at which noslip BCs are set
321 : : std::vector< std::size_t > m_noslipbcnodes;
322 : : //! Streamable surface integral nodes and normals * dA on surfaces
323 : : std::map< int, std::pair< std::vector< std::size_t >,
324 : : std::vector< tk::real > > > m_surfint;
325 : : //! Runge-Kutta stage counter
326 : : std::size_t m_stage;
327 : : //! True in the last time step
328 : : int m_finished;
329 : : //! Runge-Kutta coefficients
330 : : std::vector< tk::real > m_rkcoef;
331 : : //! Timer
332 : : std::vector< tk::Timer > m_timer;
333 : :
334 : : //! Compute number of scalar components for gradients
335 : : std::size_t ngradcomp() const;
336 : :
337 : : //! Access bound Discretization class pointer
338 : 154973 : Discretization* Disc() const {
339 [ + - ][ + - ]: 154973 : Assert( m_disc[ thisIndex ].ckLocal() != nullptr, "ckLocal() null" );
[ - + ][ - - ]
[ - - ][ - - ]
340 [ + - ][ + - ]: 154973 : return m_disc[ thisIndex ].ckLocal();
341 : : }
342 : :
343 : : //! Prepare Dirichlet boundary condition data structures
344 : : void setupDirBC( const std::vector< std::vector< int > >& cfgmask,
345 : : const std::vector< std::vector< double > >& cfgval,
346 : : std::size_t ncomp,
347 : : std::vector< std::size_t >& mask,
348 : : std::vector< double >& val );
349 : :
350 : : //! Start computing velocity divergence
351 : : void div( const tk::Fields& u, std::size_t pos = 0 );
352 : :
353 : : //! Start computing velocity gradient
354 : : void velgrad();
355 : :
356 : : //! Start computing momentum flux
357 : : void flux();
358 : :
359 : : //! Finalize computing gradient
360 : : void fingrad( tk::Fields& grad,
361 : : std::unordered_map< std::size_t, std::vector< tk::real > >& gradc );
362 : :
363 : : //! Compute local contributions to domain edge integrals
364 : : void domint();
365 : :
366 : : //! Setup lhs matrix for pressure solve
367 : : std::tuple< tk::CSR, std::vector< tk::real >, std::vector< tk::real > >
368 : : prelhs( const std::pair< std::vector< std::size_t >,
369 : : std::vector< std::size_t > >& psup );
370 : :
371 : : //! Compute chare-boundary edges
372 : : void bndEdges();
373 : :
374 : : //! Compute local contributions to boundary normals and integrals
375 : : void bndint();
376 : :
377 : : //! Combine own and communicated portions of the boundary point normals
378 : : void bnorm();
379 : :
380 : : //! Convert integrals into streamable data structures
381 : : void streamable();
382 : :
383 : : //! Generate superedge-groups for domain-edge loops
384 : : void domsuped();
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 gradients
396 : : void grad();
397 : :
398 : : //! Compute righ-hand side vector of transport equations
399 : : void rhs();
400 : :
401 : : //! Advance systems of equations
402 : : void solve();
403 : :
404 : : //! Start next time step stage
405 : : void stage();
406 : :
407 : : //! Optionally refine/derefine mesh
408 : : void refine();
409 : :
410 : : //! Compute time step size
411 : : void dt();
412 : :
413 : : //! Evaluate whether to save checkpoint/restart
414 : : void evalRestart();
415 : :
416 : : //! Apply scalar source to solution
417 : : void src();
418 : : };
419 : :
420 : : } // inciter::
|
