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-2025 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 646 : 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 : //! Initiate transfer of transfer flags (if coupled)
84 : void transferFL();
85 :
86 : //! Continue after transfer of initial conditions (if coupled)
87 : void transferIC();
88 :
89 : //! Initialize Poisson solve
90 : void pinit();
91 :
92 : //! Solve Poisson equation
93 : void psolve();
94 :
95 : //! Continue after Poisson solve
96 : void psolved();
97 :
98 : // Start time stepping
99 : void start();
100 :
101 : //! Advance equations to next time step
102 : void advance( tk::real newdt );
103 :
104 : //! Evaluate diagnostics
105 : void diag();
106 :
107 : //! Start (re-)computing domain and boundary integrals
108 : void feop();
109 :
110 : //! Receive contributions to boundary point normals on chare-boundaries
111 : void comnorm( const std::unordered_map< int,
112 : std::unordered_map< std::size_t, std::array<tk::real,4> > >& inbnd );
113 :
114 : //! Receive contributions to velocity gradients
115 : void comvgrad( const std::unordered_map< std::size_t,
116 : std::vector< tk::real > >& ingrad );
117 :
118 : //! Receive contributions to momentum flux on chare-boundaries
119 : void comflux( const std::unordered_map< std::size_t,
120 : std::vector< tk::real > >& influx );
121 :
122 : //! Receive contributions to conjugate gradients solution gradient
123 : void comsgrad( const std::unordered_map< std::size_t,
124 : std::vector< tk::real > >& ingrad );
125 :
126 : //! Receive contributions to gradient
127 : void comgrad( const std::unordered_map< std::size_t,
128 : std::vector< tk::real > >& ingrad );
129 :
130 : //! Receive contributions to right-hand side vector on chare-boundaries
131 : void comrhs( const std::unordered_map< std::size_t,
132 : std::vector< tk::real > >& inrhs );
133 :
134 : //! Receive contributions to velocity divergence on chare-boundaries
135 : void comdiv( const std::unordered_map< std::size_t, tk::real >& indiv );
136 :
137 : //! Solution has been updated
138 : void solved();
139 :
140 : //! Evaluate residuals
141 : void evalres( const std::vector< tk::real >& l2res );
142 :
143 : //! Receive new mesh from Refiner
144 : void resizePostAMR(
145 : const std::vector< std::size_t >& ginpoel,
146 : const tk::UnsMesh::Chunk& chunk,
147 : const tk::UnsMesh::Coords& coord,
148 : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
149 : const std::unordered_map< std::size_t, std::size_t >& addedTets,
150 : const std::set< std::size_t >& removedNodes,
151 : const std::unordered_map< int, std::unordered_set< std::size_t > >&
152 : nodeCommMap,
153 : const std::map< int, std::vector< std::size_t > >& bface,
154 : const std::map< int, std::vector< std::size_t > >& bnode,
155 : const std::vector< std::size_t >& triinpoel );
156 :
157 : //! Const-ref access to current solution
158 : //! \return Const-ref to current solution
159 : const tk::Fields& solution() const { return m_u; }
160 :
161 : //! Compute integral quantities for output
162 : void integrals();
163 :
164 : //! Compute recent conjugate gradients solution gradient
165 : void sgrad();
166 :
167 : //! Evaluate whether to continue with next time step
168 : void step();
169 :
170 : // Evaluate whether to do load balancing
171 : void evalLB( int nrestart );
172 :
173 : /** @name Charm++ pack/unpack serializer member functions */
174 : ///@{
175 : //! \brief Pack/Unpack serialize member function
176 : //! \param[in,out] p Charm++'s PUP::er serializer object reference
177 2026 : void pup( PUP::er &p ) override {
178 : p | m_disc;
179 : p | m_cgpre;
180 2026 : p | m_nrhs;
181 2026 : p | m_nnorm;
182 2026 : p | m_ngrad;
183 2026 : p | m_nsgrad;
184 2026 : p | m_nvgrad;
185 2026 : p | m_nflux;
186 2026 : p | m_ndiv;
187 2026 : p | m_nbpint;
188 2026 : p | m_np;
189 2026 : p | m_bnode;
190 2026 : p | m_bface;
191 2026 : p | m_triinpoel;
192 2026 : p | m_u;
193 2026 : p | m_un;
194 2026 : p | m_grad;
195 2026 : p | m_gradc;
196 2026 : p | m_vgrad;
197 2026 : p | m_vgradc;
198 2026 : p | m_flux;
199 2026 : p | m_fluxc;
200 2026 : p | m_div;
201 2026 : p | m_divc;
202 2026 : p | m_sgrad;
203 2026 : p | m_sgradc;
204 2026 : p | m_rhs;
205 2026 : p | m_rhsc;
206 : p | m_diag;
207 2026 : p | m_bnorm;
208 2026 : p | m_bnormc;
209 2026 : p | m_bndpoinint;
210 2026 : p | m_domedgeint;
211 2026 : p | m_bpint;
212 : p | m_dsupedge;
213 : p | m_dsupint;
214 2026 : p | m_dirbcmask;
215 2026 : p | m_dirbcval;
216 2026 : p | m_dirbcmaskp;
217 2026 : p | m_dirbcvalp;
218 2026 : p | m_symbcnodes;
219 2026 : p | m_symbcnorms;
220 2026 : p | m_noslipbcnodes;
221 2026 : p | m_surfint;
222 2026 : p | m_stage;
223 2026 : p | m_finished;
224 2026 : p | m_rkcoef;
225 2026 : p | m_timer;
226 2026 : }
227 : //! \brief Pack/Unpack serialize operator|
228 : //! \param[in,out] p Charm++'s PUP::er serializer object reference
229 : //! \param[in,out] i LohCG object reference
230 : friend void operator|( PUP::er& p, LohCG& i ) { i.pup(p); }
231 : //@}
232 :
233 : private:
234 : //! Discretization proxy
235 : CProxy_Discretization m_disc;
236 : //! Conjugate Gradients Charm++ proxy for pressure solve
237 : tk::CProxy_ConjugateGradients m_cgpre;
238 : //! Counter for right-hand side vector nodes updated
239 : std::size_t m_nrhs;
240 : //! Counter for receiving boundary point normals
241 : std::size_t m_nnorm;
242 : //! Counter for receiving gradient
243 : std::size_t m_ngrad;
244 : //! Counter for receiving conjugrate gradient solution gradient
245 : std::size_t m_nsgrad;
246 : //! Counter for receiving velocity gradient
247 : std::size_t m_nvgrad;
248 : //! Counter for receiving momentum flux
249 : std::size_t m_nflux;
250 : //! Counter for receiving boundary velocity divergences
251 : std::size_t m_ndiv;
252 : //! Counter for receiving boundary point integrals
253 : std::size_t m_nbpint;
254 : //! Count number of Poisson solves during setup
255 : std::size_t m_np;
256 : //! Boundary node lists mapped to side set ids used in the input file
257 : std::map< int, std::vector< std::size_t > > m_bnode;
258 : //! Boundary face lists mapped to side set ids used in the input file
259 : std::map< int, std::vector< std::size_t > > m_bface;
260 : //! Boundary triangle face connecitivity where BCs are set by user
261 : std::vector< std::size_t > m_triinpoel;
262 : //! Unknown/solution vector at mesh nodes
263 : tk::Fields m_u;
264 : //! Unknown/solution vector at mesh nodes at previous time step
265 : tk::Fields m_un;
266 : //! Gradient in mesh nodes
267 : tk::Fields m_grad;
268 : //! Gradient receive buffer
269 : std::unordered_map< std::size_t, std::vector< tk::real > > m_gradc;
270 : //! Velocity gradient in mesh nodes
271 : tk::Fields m_vgrad;
272 : //! Velocity gradient receive buffer
273 : std::unordered_map< std::size_t, std::vector< tk::real > > m_vgradc;
274 : //! Momentum flux in mesh nodes
275 : tk::Fields m_flux;
276 : //! Momentum flux receive buffer
277 : std::unordered_map< std::size_t, std::vector< tk::real > > m_fluxc;
278 : //! Velocity divergence
279 : std::vector< tk::real > m_div;
280 : //! Receive buffer for communication of the velocity divergence
281 : //! \details Key: global node id, value: velocity divergence
282 : std::unordered_map< std::size_t, tk::real > m_divc;
283 : //! Conjugate gradient solution gradient in mesh nodes
284 : tk::Fields m_sgrad;
285 : //! Conjugate gradient solution gradient receive buffer
286 : std::unordered_map< std::size_t, std::vector< tk::real > > m_sgradc;
287 : //! Right-hand side vector (for the high order system)
288 : tk::Fields m_rhs;
289 : //! Receive buffer for communication of the right hand side
290 : //! \details Key: global node id, value: rhs for all scalar components per
291 : //! node.
292 : std::unordered_map< std::size_t, std::vector< tk::real > > m_rhsc;
293 : //! Diagnostics object
294 : NodeDiagnostics m_diag;
295 : //! Boundary point normals
296 : //! \details Outer key: side set id. Inner key: global node id of boundary
297 : //! point, value: weighted normal vector, inverse distance square.
298 : std::unordered_map< int,
299 : std::unordered_map< std::size_t, std::array< tk::real, 4 > > > m_bnorm;
300 : //! Boundary point normals receive buffer
301 : //! \details Outer key: side set id. Inner key: global node id of boundary
302 : //! point, value: weighted normals and inverse distance square.
303 : decltype(m_bnorm) m_bnormc;
304 : //! Boundary point integrals
305 : //! \details Key: global node id of boundary point, value: boundary point
306 : //! integral contributions.
307 : std::unordered_map< std::size_t, std::array<tk::real,3> > m_bndpoinint;
308 : //! Domain edge integrals
309 : std::unordered_map< tk::UnsMesh::Edge, std::array< tk::real, 4 >,
310 : tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> > m_domedgeint;
311 : //! Streamable boundary point integrals
312 : std::vector< tk::real > m_bpint;
313 : //! Superedge (tet, face, edge) end points with local ids for domain edges
314 : std::array< std::vector< std::size_t >, 3 > m_dsupedge;
315 : //! Superedge (tet, face, edge) domain edge integrals
316 : std::array< std::vector< tk::real >, 3 > m_dsupint;
317 : //! Nodes and their Dirichlet BC masks
318 : std::vector< std::size_t > m_dirbcmask;
319 : //! Nodes and their Dirichlet BC values
320 : std::vector< double > m_dirbcval;
321 : //! Nodes and their pressure Dirichlet BC masks
322 : std::vector< std::size_t > m_dirbcmaskp;
323 : //! Nodes and their pressure Dirichlet BC values
324 : std::vector< double > m_dirbcvalp;
325 : //! Streamable nodes at which symmetry BCs are set
326 : std::vector< std::size_t > m_symbcnodes;
327 : //! Streamable normals at nodes at which symmetry BCs are set
328 : std::vector< tk::real > m_symbcnorms;
329 : //! Streamable nodes at which noslip BCs are set
330 : std::vector< std::size_t > m_noslipbcnodes;
331 : //! Streamable surface integral nodes and normals * dA on surfaces
332 : std::map< int, std::pair< std::vector< std::size_t >,
333 : std::vector< tk::real > > > m_surfint;
334 : //! Runge-Kutta stage counter
335 : std::size_t m_stage;
336 : //! True in the last time step
337 : int m_finished;
338 : //! Runge-Kutta coefficients
339 : std::vector< tk::real > m_rkcoef;
340 : //! Timer
341 : std::vector< tk::Timer > m_timer;
342 :
343 : //! Compute number of scalar components for gradients
344 : std::size_t ngradcomp() const;
345 :
346 : //! Access bound Discretization class pointer
347 235242 : Discretization* Disc() const {
348 : Assert( m_disc[ thisIndex ].ckLocal() != nullptr, "ckLocal() null" );
349 470484 : return m_disc[ thisIndex ].ckLocal();
350 : }
351 :
352 : //! Prepare Dirichlet boundary condition data structures
353 : void setupDirBC( const std::vector< std::vector< int > >& cfgmask,
354 : const std::vector< std::vector< double > >& cfgval,
355 : std::size_t ncomp,
356 : std::vector< std::size_t >& mask,
357 : std::vector< double >& val );
358 :
359 : //! Start computing velocity divergence
360 : void div( const tk::Fields& u, std::size_t pos = 0 );
361 :
362 : //! Start computing velocity gradient
363 : void velgrad();
364 :
365 : //! Start computing momentum flux
366 : void flux();
367 :
368 : //! Finalize computing gradient
369 : void fingrad( tk::Fields& grad,
370 : std::unordered_map< std::size_t, std::vector< tk::real > >& gradc );
371 :
372 : //! Compute local contributions to domain edge integrals
373 : void domint();
374 :
375 : //! Setup lhs matrix for pressure solve
376 : std::tuple< tk::CSR, std::vector< tk::real >, std::vector< tk::real > >
377 : prelhs( const std::pair< std::vector< std::size_t >,
378 : std::vector< std::size_t > >& psup );
379 :
380 : //! Set solution in holes (if coupled)
381 : void holeset();
382 :
383 : //! Compute chare-boundary edges
384 : void bndEdges();
385 :
386 : //! Compute local contributions to boundary normals and integrals
387 : void bndint();
388 :
389 : //! Combine own and communicated portions of the boundary point normals
390 : void bnorm();
391 :
392 : //! Prepare surface integral data strurctures
393 : void prep_surfint();
394 :
395 : //! Prepare symmetry boundary condition data structures
396 : void prep_symbc();
397 :
398 : //! Prepare no-slip boundary condition data structures
399 : void prep_noslipbc();
400 :
401 : //! Convert integrals into streamable data structures
402 : void streamable();
403 :
404 : //! Generate superedge-groups for domain-edge loops
405 : void domsuped();
406 :
407 : //! Output mesh and particle fields to files
408 : void out();
409 :
410 : //! Output mesh-based fields to file
411 : void writeFields( CkCallback cb );
412 :
413 : //! Combine own and communicated portions of the integrals
414 : void merge();
415 :
416 : //! Compute gradients
417 : void grad();
418 :
419 : //! Compute righ-hand side vector of transport equations
420 : void rhs();
421 :
422 : //! Advance systems of equations
423 : void solve();
424 :
425 : //! Start next time step stage
426 : void stage();
427 :
428 : //! Optionally refine/derefine mesh
429 : void refine();
430 :
431 : //! Compute time step size
432 : void dt();
433 :
434 : //! Evaluate whether to save checkpoint/restart
435 : void evalRestart();
436 :
437 : //! Apply scalar source to solution
438 : void src();
439 : };
440 :
441 : } // inciter::
|
