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