Branch data 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-2024 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 : 772 : 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 : : //! Evaluate diagnostics
107 : : void diag();
108 : :
109 : : //! Start (re-)computing domain and boundary integrals
110 : : void feop();
111 : :
112 : : //! Receive contributions to boundary point normals on chare-boundaries
113 : : void comnorm( const std::unordered_map< int,
114 : : std::unordered_map< std::size_t, std::array<tk::real,4> > >& inbnd );
115 : :
116 : : //! Receive contributions to velocity gradients
117 : : void comvgrad( const std::unordered_map< std::size_t,
118 : : std::vector< tk::real > >& ingrad );
119 : :
120 : : //! Receive contributions to momentum flux on chare-boundaries
121 : : void comflux( const std::unordered_map< std::size_t,
122 : : std::vector< tk::real > >& influx );
123 : :
124 : : //! Receive contributions to conjugate gradients solution gradient
125 : : void comsgrad( const std::unordered_map< std::size_t,
126 : : std::vector< tk::real > >& ingrad );
127 : :
128 : : //! Receive contributions to pressure gradient
129 : : void compgrad( const std::unordered_map< std::size_t,
130 : : std::vector< tk::real > >& ingrad );
131 : :
132 : : //! Receive contributions to right-hand side vector on chare-boundaries
133 : : void comrhs( const std::unordered_map< std::size_t,
134 : : std::vector< tk::real > >& inrhs );
135 : :
136 : : //! Receive contributions to velocity divergence on chare-boundaries
137 : : void comdiv( const std::unordered_map< std::size_t, tk::real >& indiv );
138 : :
139 : : //! Receive antidiffusive and low-order contributions on chare-boundaries
140 : : void comaec( const std::unordered_map< std::size_t,
141 : : std::vector< tk::real > >& inaec );
142 : :
143 : :
144 : : //! Receive allowed limits contributions on chare-boundaries
145 : : void comalw( const std::unordered_map< std::size_t,
146 : : std::vector< tk::real > >& inalw );
147 : :
148 : : //! Receive limited antidiffusive contributions on chare-boundaries
149 : : void comlim( const std::unordered_map< std::size_t,
150 : : std::vector< tk::real > >& inlim );
151 : :
152 : : //! Evaluate residuals
153 : : void evalres( const std::vector< tk::real >& l2res );
154 : :
155 : : //! Receive new mesh from Refiner
156 : : void resizePostAMR(
157 : : const std::vector< std::size_t >& ginpoel,
158 : : const tk::UnsMesh::Chunk& chunk,
159 : : const tk::UnsMesh::Coords& coord,
160 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
161 : : const std::unordered_map< std::size_t, std::size_t >& addedTets,
162 : : const std::set< std::size_t >& removedNodes,
163 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
164 : : nodeCommMap,
165 : : const std::map< int, std::vector< std::size_t > >& bface,
166 : : const std::map< int, std::vector< std::size_t > >& bnode,
167 : : const std::vector< std::size_t >& triinpoel );
168 : :
169 : : //! Const-ref access to current solution
170 : : //! \return Const-ref to current solution
171 : : const tk::Fields& solution() const { return m_u; }
172 : :
173 : : //! Compute integral quantities for output
174 : : void integrals();
175 : :
176 : : //! Compute recent conjugate gradients solution gradient
177 : : void sgrad();
178 : :
179 : : //! Evaluate whether to continue with next time step
180 : : void step();
181 : :
182 : : // Evaluate whether to do load balancing
183 : : void evalLB( int nrestart );
184 : :
185 : : /** @name Charm++ pack/unpack serializer member functions */
186 : : ///@{
187 : : //! \brief Pack/Unpack serialize member function
188 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
189 : 2384 : void pup( PUP::er &p ) override {
190 : : p | m_disc;
191 : : p | m_cgpre;
192 : : p | m_cgmom;
193 : 2384 : p | m_nrhs;
194 : 2384 : p | m_nnorm;
195 : 2384 : p | m_naec;
196 : 2384 : p | m_nalw;
197 : 2384 : p | m_nlim;
198 : 2384 : p | m_nsgrad;
199 : 2384 : p | m_npgrad;
200 : 2384 : p | m_nvgrad;
201 : 2384 : p | m_nflux;
202 : 2384 : p | m_ndiv;
203 : 2384 : p | m_nbpint;
204 : 2384 : p | m_np;
205 : 2384 : p | m_bnode;
206 : 2384 : p | m_bface;
207 : 2384 : p | m_triinpoel;
208 : 2384 : p | m_u;
209 : 2384 : p | m_un;
210 : 2384 : p | m_pr;
211 : 2384 : p | m_p;
212 : 2384 : p | m_pc;
213 : 2384 : p | m_q;
214 : 2384 : p | m_qc;
215 : 2384 : p | m_a;
216 : 2384 : p | m_ac;
217 : 2384 : p | m_rhs;
218 : 2384 : p | m_rhsc;
219 : 2384 : p | m_sgrad;
220 : 2384 : p | m_sgradc;
221 : 2384 : p | m_pgrad;
222 : 2384 : p | m_pgradc;
223 : 2384 : p | m_vgrad;
224 : 2384 : p | m_vgradc;
225 : 2384 : p | m_flux;
226 : 2384 : p | m_fluxc;
227 : 2384 : p | m_div;
228 : 2384 : p | m_divc;
229 : : p | m_diag;
230 : 2384 : p | m_bnorm;
231 : 2384 : p | m_bnormc;
232 : 2384 : p | m_bndpoinint;
233 : 2384 : p | m_domedgeint;
234 : 2384 : p | m_bpint;
235 : : p | m_dsupedge;
236 : : p | m_dsupint;
237 : : p | m_dsuplim;
238 : 2384 : p | m_dirbcmask;
239 : 2384 : p | m_dirbcval;
240 : 2384 : p | m_dirbcmaskp;
241 : 2384 : p | m_dirbcvalp;
242 : 2384 : p | m_symbcnodes;
243 : 2384 : p | m_symbcnorms;
244 : 2384 : p | m_noslipbcnodes;
245 : 2384 : p | m_surfint;
246 : 2384 : p | m_stage;
247 : 2384 : p | m_finished;
248 : 2384 : p | m_rkcoef;
249 : 2384 : p | m_timer;
250 : 2384 : }
251 : : //! \brief Pack/Unpack serialize operator|
252 : : //! \param[in,out] p Charm++'s PUP::er serializer object reference
253 : : //! \param[in,out] i ChoCG object reference
254 : : friend void operator|( PUP::er& p, ChoCG& i ) { i.pup(p); }
255 : : //@}
256 : :
257 : : private:
258 : : //! Discretization proxy
259 : : CProxy_Discretization m_disc;
260 : : //! Conjugate Gradients Charm++ proxy for pressure solve
261 : : tk::CProxy_ConjugateGradients m_cgpre;
262 : : //! Conjugate Gradients Charm++ proxy for momentum solve
263 : : tk::CProxy_ConjugateGradients m_cgmom;
264 : : //! Counter for right-hand side vector nodes updated
265 : : std::size_t m_nrhs;
266 : : //! Counter for receiving boundary point normals
267 : : std::size_t m_nnorm;
268 : : //! Counter for receiving antidiffusive contributions
269 : : std::size_t m_naec;
270 : : //! Counter for receiving allowed limits
271 : : std::size_t m_nalw;
272 : : //! Counter for receiving limited antidiffusive contributions
273 : : std::size_t m_nlim;
274 : : //! Counter for receiving conjugrate gradient solution gradient
275 : : std::size_t m_nsgrad;
276 : : //! Counter for receiving pressure gradient
277 : : std::size_t m_npgrad;
278 : : //! Counter for receiving velocity gradient
279 : : std::size_t m_nvgrad;
280 : : //! Counter for receiving momentum flux
281 : : std::size_t m_nflux;
282 : : //! Counter for receiving boundary velocity divergences
283 : : std::size_t m_ndiv;
284 : : //! Counter for receiving boundary point integrals
285 : : std::size_t m_nbpint;
286 : : //! Count number of Poisson solves during setup
287 : : std::size_t m_np;
288 : : //! Boundary node lists mapped to side set ids used in the input file
289 : : std::map< int, std::vector< std::size_t > > m_bnode;
290 : : //! Boundary face lists mapped to side set ids used in the input file
291 : : std::map< int, std::vector< std::size_t > > m_bface;
292 : : //! Boundary triangle face connecitivity where BCs are set by user
293 : : std::vector< std::size_t > m_triinpoel;
294 : : //! Unknown/solution vector at mesh nodes
295 : : tk::Fields m_u;
296 : : //! Unknown/solution vector at mesh nodes at previous time step
297 : : tk::Fields m_un;
298 : : //! Pressure
299 : : std::vector< tk::real > m_pr;
300 : : //! Max/min antidiffusive edge contributions at mesh nodes
301 : : tk::Fields m_p;
302 : : //! Receive buffer for max/min antidiffusive edge contributions
303 : : //! \details Key: global node id, value: max/min antidiff edge contributions
304 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_pc;
305 : : //! Max/min allowed limits at mesh nodes
306 : : tk::Fields m_q;
307 : : //! Receive buffer for max/min allowed limits
308 : : //! \details Key: global node id, value: max/min allowed limits in nodes.
309 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_qc;
310 : : //! Limited antidiffusive contributions at mesh nodes
311 : : tk::Fields m_a;
312 : : //! Receive buffer for limited antidiffusive contributions
313 : : //! \details Key: global node id, value: limited antidiffusive contributions
314 : : //! in nodes.
315 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_ac;
316 : : //! Right-hand side vector (for the high order system)
317 : : tk::Fields m_rhs;
318 : : //! Receive buffer for communication of the right hand side
319 : : //! \details Key: global node id, value: rhs for all scalar components per
320 : : //! node.
321 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_rhsc;
322 : : //! Conjugate gradient solution gradient in mesh nodes
323 : : tk::Fields m_sgrad;
324 : : //! Conjugate gradient solution gradient receive buffer
325 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_sgradc;
326 : : //! Pressure gradient in mesh nodes
327 : : tk::Fields m_pgrad;
328 : : //! Pressure gradient receive buffer
329 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_pgradc;
330 : : //! Velocity gradient in mesh nodes
331 : : tk::Fields m_vgrad;
332 : : //! Velocity gradient receive buffer
333 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_vgradc;
334 : : //! Momentum flux in mesh nodes
335 : : tk::Fields m_flux;
336 : : //! Momentum flux receive buffer
337 : : std::unordered_map< std::size_t, std::vector< tk::real > > m_fluxc;
338 : : //! Velocity divergence
339 : : std::vector< tk::real > m_div;
340 : : //! Receive buffer for communication of the velocity divergence
341 : : //! \details Key: global node id, value: velocity divergence
342 : : std::unordered_map< std::size_t, tk::real > m_divc;
343 : : //! Diagnostics object
344 : : NodeDiagnostics m_diag;
345 : : //! Boundary point normals
346 : : //! \details Outer key: side set id. Inner key: global node id of boundary
347 : : //! point, value: weighted normal vector, inverse distance square.
348 : : std::unordered_map< int,
349 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > > m_bnorm;
350 : : //! Boundary point normals receive buffer
351 : : //! \details Outer key: side set id. Inner key: global node id of boundary
352 : : //! point, value: weighted normals and inverse distance square.
353 : : decltype(m_bnorm) m_bnormc;
354 : : //! Boundary point integrals
355 : : //! \details Key: global node id of boundary point, value: boundary point
356 : : //! integral contributions.
357 : : std::unordered_map< std::size_t, std::array<tk::real,3> > m_bndpoinint;
358 : : //! Domain edge integrals
359 : : std::unordered_map< tk::UnsMesh::Edge, std::array< tk::real, 5 >,
360 : : tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> > m_domedgeint;
361 : : //! Streamable boundary point integrals
362 : : std::vector< tk::real > m_bpint;
363 : : //! Superedge (tet, face, edge) end points with local ids for domain edges
364 : : std::array< std::vector< std::size_t >, 3 > m_dsupedge;
365 : : //! Superedge (tet, face, edge) domain edge integrals
366 : : std::array< std::vector< tk::real >, 3 > m_dsupint;
367 : : //! FCT limiter coefficients in domain superedges
368 : : std::array< std::vector< tk::real >, 3 > m_dsuplim;
369 : : //! Nodes and their Dirichlet BC masks
370 : : std::vector< std::size_t > m_dirbcmask;
371 : : //! Nodes and their Dirichlet BC values
372 : : std::vector< double > m_dirbcval;
373 : : //! Nodes and their pressure Dirichlet BC masks
374 : : std::vector< std::size_t > m_dirbcmaskp;
375 : : //! Nodes and their pressure Dirichlet BC values
376 : : std::vector< double > m_dirbcvalp;
377 : : //! Streamable nodes at which symmetry BCs are set
378 : : std::vector< std::size_t > m_symbcnodes;
379 : : //! Streamable normals at nodes at which symmetry BCs are set
380 : : std::vector< tk::real > m_symbcnorms;
381 : : //! Streamable nodes at which noslip BCs are set
382 : : std::vector< std::size_t > m_noslipbcnodes;
383 : : //! Streamable surface integral nodes and normals * dA on surfaces
384 : : std::map< int, std::pair< std::vector< std::size_t >,
385 : : std::vector< tk::real > > > m_surfint;
386 : : //! Runge-Kutta stage counter
387 : : std::size_t m_stage;
388 : : //! True in the last time step
389 : : int m_finished;
390 : : //! Runge-Kutta coefficients
391 : : std::vector< tk::real > m_rkcoef;
392 : : //! Timer
393 : : std::vector< tk::Timer > m_timer;
394 : :
395 : : //! Access bound Discretization class pointer
396 : 427276 : Discretization* Disc() const {
397 : : Assert( m_disc[ thisIndex ].ckLocal() != nullptr, "ckLocal() null" );
398 : 854552 : return m_disc[ thisIndex ].ckLocal();
399 : : }
400 : :
401 : : //! Access bound momentum solve matrix
402 : 40 : tk::CSR& Lhs() {
403 : : Assert( m_cgmom[ thisIndex ].ckLocal() != nullptr, "ckLocal() null" );
404 : 80 : return m_cgmom[ thisIndex ].ckLocal()->lhs();
405 : : }
406 : :
407 : : //! Prepare Dirichlet boundary condition data structures
408 : : void setupDirBC( const std::vector< std::vector< int > >& cfgmask,
409 : : const std::vector< std::vector< double > >& cfgval,
410 : : std::size_t ncomp,
411 : : std::vector< std::size_t >& mask,
412 : : std::vector< double >& val );
413 : :
414 : : //! Start computing velocity divergence
415 : : void div( const tk::Fields& u );
416 : :
417 : : //! Start computing velocity gradient
418 : : void velgrad();
419 : :
420 : : //! Start computing momentum flux
421 : : void flux();
422 : :
423 : : //! Finalize computing gradient
424 : : void fingrad( tk::Fields& grad,
425 : : std::unordered_map< std::size_t, std::vector< tk::real > >& gradc );
426 : :
427 : : //! Finalize computing pressure gradient
428 : : void finpgrad();
429 : :
430 : : //! Compute pressure gradient
431 : : void pgrad();
432 : :
433 : : //! Compute local contributions to domain edge integrals
434 : : void domint();
435 : :
436 : : //! Setup lhs matrix for pressure solve
437 : : std::tuple< tk::CSR, std::vector< tk::real >, std::vector< tk::real > >
438 : : prelhs( const std::pair< std::vector< std::size_t >,
439 : : std::vector< std::size_t > >& psup );
440 : :
441 : : //! Setup empty lhs matrix for momentum solve
442 : : std::tuple< tk::CSR, std::vector< tk::real >, std::vector< tk::real > >
443 : : momlhs( const std::pair< std::vector< std::size_t >,
444 : : std::vector< std::size_t > >& psup );
445 : :
446 : : //! Compute chare-boundary edges
447 : : void bndEdges();
448 : :
449 : : //! Compute local contributions to boundary normals and integrals
450 : : void bndint();
451 : :
452 : : //! Combine own and communicated portions of the boundary point normals
453 : : void bnorm();
454 : :
455 : : //! Convert integrals into streamable data structures
456 : : void streamable();
457 : :
458 : : //! Generate superedge-groups for domain-edge loops
459 : : void domsuped();
460 : :
461 : : //! Output mesh and particle fields to files
462 : : void out();
463 : :
464 : : //! Output mesh-based fields to file
465 : : void writeFields( CkCallback cb );
466 : :
467 : : //! Combine own and communicated portions of the integrals
468 : : void merge();
469 : :
470 : : //! Fill lhs matrix of transport equations
471 : : void lhs();
472 : :
473 : : //! Compute righ-hand side vector of transport equations
474 : : void rhs();
475 : :
476 : : //! Continue with flux-corrected transport if enabled
477 : : void fct();
478 : :
479 : : //! Compute antidiffusive contributions: P+/-
480 : : void aec();
481 : :
482 : : //! Compute allowed limits, Q+/-
483 : : void alw();
484 : :
485 : : //! Compute limit coefficients
486 : : void lim();
487 : :
488 : : //! Advance systems of equations
489 : : void solve();
490 : :
491 : : //! Compute advective-diffusive prediction of momentum/transport
492 : : void pred();
493 : :
494 : : //! Compute pressure correction
495 : : void corr();
496 : :
497 : : //! Optionally refine/derefine mesh
498 : : void refine();
499 : :
500 : : //! Compute time step size
501 : : void dt();
502 : :
503 : : //! Evaluate whether to save checkpoint/restart
504 : : void evalRestart();
505 : :
506 : : //! Apply boundary conditions
507 : : void BC( tk::Fields& u, tk::real t );
508 : :
509 : : //! Apply scalar source to solution
510 : : void src();
511 : : };
512 : :
513 : : } // inciter::
|
