Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/DistFCT.cpp
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 Charm++ chare array for distributed flux-corrected transport
9 : : \details Charm++ chare array for asynchronous distributed
10 : : flux-corrected transport (FCT).
11 : : */
12 : : // *****************************************************************************
13 : :
14 : : #include <string>
15 : : #include <cmath>
16 : : #include <array>
17 : : #include <set>
18 : : #include <algorithm>
19 : :
20 : : #include "ContainerUtil.hpp"
21 : : #include "DistFCT.hpp"
22 : :
23 : : namespace inciter {
24 : :
25 : : extern ctr::InputDeck g_inputdeck;
26 : :
27 : : } // inciter::
28 : :
29 : : using inciter::DistFCT;
30 : :
31 : 654 : DistFCT::DistFCT( int nchare,
32 : : std::size_t nu,
33 : : std::size_t np,
34 : : const tk::NodeCommMap& nodeCommMap,
35 : : const std::unordered_map< std::size_t, std::size_t >& bid,
36 : : const std::unordered_map< std::size_t, std::size_t >& lid,
37 : 654 : const std::vector< std::size_t >& inpoel ) :
38 : : m_naec( 0 ),
39 : : m_nalw( 0 ),
40 : : m_nlim( 0 ),
41 : 654 : m_nchare( static_cast< std::size_t >( nchare ) ),
42 : : m_nodeCommMap( nodeCommMap ),
43 : : m_bid( bid ),
44 : : m_lid( lid ),
45 : : m_inpoel( inpoel ),
46 : : m_fluxcorrector( m_inpoel.size() ),
47 : : m_p( nu, np*2 ),
48 : : m_q( nu, np*2 ),
49 : : m_a( nu, np ),
50 : : m_pc(),
51 : : m_qc(),
52 : : m_ac(),
53 : : m_ul(),
54 : : m_dul(),
55 [ + - ][ + - ]: 654 : m_du()
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
56 : : // *****************************************************************************
57 : : // Constructor
58 : : //! \param[in] nchare Total number of worker chares
59 : : //! \param[in] nu Number of unknowns in solution vector
60 : : //! \param[in] np Total number of properties, i.e., scalar variables or
61 : : //! components, per unknown in solution vector
62 : : //! \param[in] nodeCommMap Global mesh node IDs associated to chare IDs
63 : : //! bordering the mesh chunk we operate on
64 : : //! \param[in] bid Local chare-boundary mesh node IDs at which we receive
65 : : //! contributions associated to global mesh node IDs of mesh elements we
66 : : //! contribute to
67 : : //! \param[in] lid Local mesh node ids associated to the global ones of owned
68 : : //! elements
69 : : //! \param[in] inpoel Mesh connectivity of our chunk of the mesh
70 : : // *****************************************************************************
71 : : {
72 [ + - ]: 654 : resizeComm(); // Size communication buffers
73 : 654 : }
74 : :
75 : : void
76 : 730 : DistFCT::resizeComm()
77 : : // *****************************************************************************
78 : : // Size FCT communication buffers
79 : : //! \details The size of the communication buffers are determined based on
80 : : //! m_bid.size() and m_a.nprop().
81 : : // *****************************************************************************
82 : : {
83 : 730 : auto bs = m_bid.size();
84 : 730 : auto np = m_a.nprop();
85 : :
86 : 730 : m_pc.resize( bs );
87 [ + + ][ + - ]: 91650 : for (auto& b : m_pc) b.resize( np*2 );
88 : 730 : m_qc.resize( bs );
89 [ + + ][ + - ]: 91650 : for (auto& b : m_qc) b.resize( np*2 );
90 : 730 : m_ac.resize( bs );
91 [ + + ][ + - ]: 91650 : for (auto& b : m_ac) b.resize( np );
92 : 730 : }
93 : :
94 : : void
95 : 47 : DistFCT::remap( const Discretization& d )
96 : : // *****************************************************************************
97 : : // Remap local ids after a mesh node reorder
98 : : //! \param[in] d Discretization proxy to read mesh data from
99 : : // *****************************************************************************
100 : : {
101 : 47 : m_lid = d.Lid();
102 : 47 : m_inpoel = d.Inpoel();
103 : 47 : }
104 : :
105 : : void
106 : 76 : DistFCT::resize( std::size_t nu,
107 : : const tk::NodeCommMap& nodeCommMap,
108 : : const std::unordered_map< std::size_t, std::size_t >& bid,
109 : : const std::unordered_map< std::size_t, std::size_t >& lid,
110 : : const std::vector< std::size_t >& inpoel )
111 : : // *****************************************************************************
112 : : // Resize FCT data structures (e.g., after mesh refinement)
113 : : //! \param[in] nu New number of unknowns in solution vector
114 : : //! \param[in] nodeCommMap New global mesh node IDs associated to chare IDs
115 : : //! bordering the mesh chunk we operate on
116 : : //! \param[in] bid New local chare-boundary mesh node IDs at which we receive
117 : : //! contributions associated to global mesh node IDs of mesh elements we
118 : : //! contribute to
119 : : //! \param[in] lid New local mesh node ids associated to the global ones of
120 : : //! owned elements
121 : : //! \param[in] inpoel Mesh connectivity of our chunk of the mesh
122 : : // *****************************************************************************
123 : : {
124 : 76 : m_nodeCommMap = nodeCommMap;
125 : 76 : m_bid = bid;
126 : 76 : m_lid = lid;
127 : 76 : m_inpoel = inpoel;
128 : :
129 : 76 : m_p.resize( nu );
130 : 76 : m_q.resize( nu );
131 : 76 : m_a.resize( nu );
132 : 76 : resizeComm();
133 : :
134 : 76 : m_fluxcorrector.resize( m_inpoel.size() );
135 : 76 : }
136 : :
137 : : tk::Fields
138 : 18283 : DistFCT::diff( const Discretization& d, const tk::Fields& Un )
139 : : // *****************************************************************************
140 : : // Compute mass diffusion rhs contribution required for the low order solution
141 : : //! \param[in] d Discretization proxy to read mesh data from
142 : : //! \param[in] Un Solution at the previous time step
143 : : //! \return Mass diffusion contribution to the RHS of the low order system
144 : : // *****************************************************************************
145 : : {
146 : 18283 : return m_fluxcorrector.diff( d.Coord(), m_inpoel, Un );
147 : : }
148 : :
149 : : void
150 : 18283 : DistFCT::next()
151 : : // *****************************************************************************
152 : : // Prepare for next time step stage
153 : : // *****************************************************************************
154 : : {
155 : : // Initialize FCT data structures for new time step
156 : 18283 : m_p.fill( 0.0 );
157 : 18283 : m_a.fill( 0.0 );
158 [ + + ]: 2755420 : for (std::size_t p=0; p<m_q.nunk(); ++p)
159 [ + + ]: 7596754 : for (ncomp_t c=0; c<m_q.nprop()/2; ++c) {
160 : 4859617 : m_q(p,c*2+0,0) = -std::numeric_limits< tk::real >::max();
161 : 4859617 : m_q(p,c*2+1,0) = std::numeric_limits< tk::real >::max();
162 : : }
163 : :
164 [ + + ][ + - ]: 1042826 : for (auto& b : m_pc) std::fill( begin(b), end(b), 0.0 );
165 [ + + ][ + - ]: 1042826 : for (auto& b : m_ac) std::fill( begin(b), end(b), 0.0 );
166 [ + + ]: 1042826 : for (auto& b : m_qc)
167 [ + + ]: 3131826 : for (ncomp_t c=0; c<m_q.nprop()/2; ++c) {
168 : 2107283 : b[c*2+0] = -std::numeric_limits< tk::real >::max();
169 : 2107283 : b[c*2+1] = std::numeric_limits< tk::real >::max();
170 : : }
171 : :
172 [ + - ]: 18283 : thisProxy[ thisIndex ].wait4fct();
173 [ + - ]: 18283 : thisProxy[ thisIndex ].wait4app();
174 : 18283 : }
175 : :
176 : : void
177 : 18283 : DistFCT::aec(
178 : : const Discretization& d,
179 : : const tk::Fields& dUh,
180 : : const tk::Fields& Un,
181 : : const std::unordered_map< std::size_t,
182 : : std::vector< std::pair< bool, tk::real > > >& bcdir,
183 : : const std::unordered_map< int,
184 : : std::unordered_set< std::size_t > >& symbcnodemap,
185 : : const std::unordered_map< int,
186 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > >& bnorm )
187 : : // *****************************************************************************
188 : : // Compute and sum antidiffusive element contributions (AEC) to mesh nodes
189 : : //! \param[in] d Discretization proxy to read mesh data from
190 : : //! \param[in] dUh Increment of the high order solution
191 : : //! \param[in] Un Solution at the previous time step
192 : : //! \param[in] bcdir Vector of pairs of bool and boundary condition value
193 : : //! associated to mesh node IDs at which to set Dirichlet boundary conditions.
194 : : //! Note that this BC data structure must include boundary conditions set
195 : : //! across all PEs, not just the ones need to be set on this PE.
196 : : //! \param[in] symbcnodemap Unique set of node ids at which to set symmetry BCs
197 : : //! \param[in] bnorm Face normals in boundary points: key global node id,
198 : : //! value: unit normal, outer key: side set id
199 : : //! \details This function computes and starts communicating m_p, which stores
200 : : //! the sum of all positive (negative) antidiffusive element contributions to
201 : : //! nodes (Lohner: P^{+,-}_i), see also FluxCorrector::aec().
202 : : // *****************************************************************************
203 : : {
204 : : // Store a copy of the high order solution increment for later
205 : 18283 : m_du = dUh;
206 : :
207 : : // Compute and sum antidiffusive element contributions to mesh nodes. Note
208 : : // that the sums are complete on nodes that are not shared with other chares
209 : : // and only partial sums on chare-boundary nodes.
210 : 18283 : m_fluxcorrector.aec( d.Coord(), m_inpoel, d.Vol(), bcdir, symbcnodemap, bnorm,
211 : 18283 : Un, m_p );
212 : :
213 [ + + ]: 18283 : if (d.NodeCommMap().empty())
214 : 499 : comaec_complete();
215 : : else // send contributions to chare-boundary nodes to fellow chares
216 [ + + ]: 183666 : for (const auto& [c,n] : d.NodeCommMap()) {
217 [ + - ]: 165882 : std::vector< std::vector< tk::real > > p( n.size() );
218 : 165882 : std::size_t j = 0;
219 [ + + ][ + - ]: 1769170 : for (auto i : n) p[ j++ ] = m_p[ tk::cref_find(m_lid,i) ];
[ + - ]
220 [ + - ][ + - ]: 165882 : thisProxy[ c ].comaec( std::vector<std::size_t>(begin(n),end(n)), p );
[ + - ]
221 : : }
222 : :
223 : 18283 : ownaec_complete( bcdir );
224 : 18283 : }
225 : :
226 : : void
227 : 165882 : DistFCT::comaec( const std::vector< std::size_t >& gid,
228 : : const std::vector< std::vector< tk::real > >& P )
229 : : // *****************************************************************************
230 : : // Receive sums of antidiffusive element contributions on chare-boundaries
231 : : //! \param[in] gid Global mesh node IDs at which we receive AEC contributions
232 : : //! \param[in] P Partial sums of positive (negative) antidiffusive element
233 : : //! contributions to chare-boundary nodes
234 : : //! \details This function receives contributions to m_p, which stores the
235 : : //! sum of all positive (negative) antidiffusive element contributions to
236 : : //! nodes (Lohner: P^{+,-}_i), see also FluxCorrector::aec(). While m_p stores
237 : : //! own contributions, m_pc collects the neighbor chare contributions during
238 : : //! communication. This way work on m_p and m_pc is overlapped. The two are
239 : : //! combined in lim().
240 : : // *****************************************************************************
241 : : {
242 [ - + ][ - - ]: 165882 : Assert( P.size() == gid.size(), "Size mismatch" );
[ - - ][ - - ]
243 : :
244 : : using tk::operator+=;
245 : :
246 [ + + ]: 1769170 : for (std::size_t i=0; i<gid.size(); ++i) {
247 : 1603288 : auto bid = tk::cref_find( m_bid, gid[i] );
248 [ - + ][ - - ]: 1603288 : Assert( bid < m_pc.size(), "Indexing out of bounds" );
[ - - ][ - - ]
249 : 1603288 : m_pc[ bid ] += P[i];
250 : : }
251 : :
252 [ + + ]: 165882 : if (++m_naec == m_nodeCommMap.size()) {
253 : 17784 : m_naec = 0;
254 : 17784 : comaec_complete();
255 : : }
256 : 165882 : }
257 : :
258 : : void
259 : 18283 : DistFCT::alw( const tk::Fields& Un,
260 : : const tk::Fields& Ul,
261 : : tk::Fields&& dUl,
262 : : const CProxy_DiagCG& host )
263 : : // *****************************************************************************
264 : : // Compute the maximum and minimum unknowns of elements surrounding nodes
265 : : //! \param[in] Un Solution at the previous time step
266 : : //! \param[in] Ul Low order solution
267 : : //! \param[in] dUl Low order solution increment
268 : : //! \param[in] host DiagCG Charm++ proxy we interoperate with
269 : : //! \details This function computes and starts communicating m_q, which stores
270 : : //! the maximum and mimimum unknowns of all elements surrounding each node
271 : : //! (Lohner: u^{max,min}_i), see also FluxCorrector::alw().
272 : : // *****************************************************************************
273 : : {
274 : : // Store a copy of the low order solution vector and its increment for later
275 : 18283 : m_ul = Ul;
276 : 18283 : m_dul = std::move(dUl);
277 : :
278 : : // Store discretization scheme proxy
279 : 18283 : m_host = host;
280 : :
281 : : // Compute the maximum and minimum unknowns of all elements surrounding nodes
282 : : // Note that the maximum and minimum unknowns are complete on nodes that are
283 : : // not shared with other chares and only partially complete on chare-boundary
284 : : // nodes.
285 : 18283 : m_fluxcorrector.alw( m_inpoel, Un, Ul, m_q );
286 : :
287 [ + + ]: 18283 : if (m_nodeCommMap.empty())
288 : 499 : comalw_complete();
289 : : else // send contributions at chare-boundary nodes to fellow chares
290 [ + + ]: 183666 : for (const auto& [c,n] : m_nodeCommMap) {
291 [ + - ]: 165882 : std::vector< std::vector< tk::real > > q( n.size() );
292 : 165882 : std::size_t j = 0;
293 [ + + ][ + - ]: 1769170 : for (auto i : n) q[ j++ ] = m_q[ tk::cref_find(m_lid,i) ];
[ + - ]
294 [ + - ][ + - ]: 165882 : thisProxy[ c ].comalw( std::vector<std::size_t>(begin(n),end(n)), q );
[ + - ]
295 : : }
296 : :
297 : 18283 : ownalw_complete();
298 : 18283 : }
299 : :
300 : : void
301 : 165882 : DistFCT::comalw( const std::vector< std::size_t >& gid,
302 : : const std::vector< std::vector< tk::real > >& Q )
303 : : // *****************************************************************************
304 : : // Receive contributions to the maxima and minima of unknowns of all elements
305 : : // surrounding mesh nodes on chare-boundaries
306 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
307 : : //! \param[in] Q Partial contributions to maximum and minimum unknowns of all
308 : : //! elements surrounding nodes to chare-boundary nodes
309 : : //! \details This function receives contributions to m_q, which stores the
310 : : //! maximum and mimimum unknowns of all elements surrounding each node
311 : : //! (Lohner: u^{max,min}_i), see also FluxCorrector::alw(). While m_q stores
312 : : //! own contributions, m_qc collects the neighbor chare contributions during
313 : : //! communication. This way work on m_q and m_qc is overlapped. The two are
314 : : //! combined in lim().
315 : : // *****************************************************************************
316 : : {
317 [ - + ][ - - ]: 165882 : Assert( Q.size() == gid.size(), "Size mismatch" );
[ - - ][ - - ]
318 : :
319 [ + + ]: 1769170 : for (std::size_t i=0; i<gid.size(); ++i) {
320 : 1603288 : auto bid = tk::cref_find( m_bid, gid[i] );
321 [ - + ][ - - ]: 1603288 : Assert( bid < m_qc.size(), "Indexing out of bounds" );
[ - - ][ - - ]
322 : 1603288 : auto& o = m_qc[ bid ];
323 : 1603288 : const auto& q = Q[i];
324 [ + + ]: 5469800 : for (ncomp_t c=0; c<m_q.nprop()/2; ++c) {
325 [ + + ]: 3866512 : if (q[c*2+0] > o[c*2+0]) o[c*2+0] = q[c*2+0];
326 [ + + ]: 3866512 : if (q[c*2+1] < o[c*2+1]) o[c*2+1] = q[c*2+1];
327 : : }
328 : : }
329 : :
330 [ + + ]: 165882 : if (++m_nalw == m_nodeCommMap.size()) {
331 : 17784 : m_nalw = 0;
332 : 17784 : comalw_complete();
333 : : }
334 : 165882 : }
335 : :
336 : : void
337 : 18283 : DistFCT::lim( const std::unordered_map< std::size_t,
338 : : std::vector< std::pair< bool, tk::real > > >& bcdir )
339 : : // *****************************************************************************
340 : : // Compute the limited antidiffusive element contributions
341 : : //! \param[in] bcdir Vector of pairs of bool and boundary condition value
342 : : //! associated to mesh node IDs at which to set Dirichlet boundary conditions.
343 : : //! Note that this BC data structure must include boundary conditions set
344 : : //! across all PEs, not just the ones need to be set on this PE.
345 : : //! \details This function computes and starts communicating m_a, which stores
346 : : //! the limited antidiffusive element contributions assembled to nodes
347 : : //! (Lohner: AEC^c), see also FluxCorrector::limit().
348 : : // *****************************************************************************
349 : : {
350 : 18283 : m_fluxcorrector.verify( m_nchare, m_inpoel, m_du, m_dul );
351 : :
352 : : // Combine own and communicated contributions to P and Q
353 [ + + ]: 1042826 : for (const auto& b : m_bid) {
354 [ + - ]: 1024543 : auto lid = tk::cref_find( m_lid, b.first );
355 : 1024543 : const auto& bpc = m_pc[ b.second ];
356 : 1024543 : const auto& bqc = m_qc[ b.second ];
357 [ + + ]: 3131826 : for (ncomp_t c=0; c<m_p.nprop()/2; ++c) {
358 [ + - ]: 2107283 : m_p(lid,c*2+0,0) += bpc[c*2+0];
359 [ + - ]: 2107283 : m_p(lid,c*2+1,0) += bpc[c*2+1];
360 [ + - ][ + + ]: 2107283 : if (bqc[c*2+0] > m_q(lid,c*2+0,0)) m_q(lid,c*2+0,0) = bqc[c*2+0];
[ + - ]
361 [ + - ][ + + ]: 2107283 : if (bqc[c*2+1] < m_q(lid,c*2+1,0)) m_q(lid,c*2+1,0) = bqc[c*2+1];
[ + - ]
362 : : }
363 : : }
364 : :
365 : 18283 : m_fluxcorrector.lim( m_inpoel, bcdir, m_p, m_ul, m_q, m_a );
366 : :
367 [ + + ]: 18283 : if (m_nodeCommMap.empty())
368 : 499 : comlim_complete();
369 : : else // send contributions to chare-boundary nodes to fellow chares
370 [ + + ]: 183666 : for (const auto& [c,n] : m_nodeCommMap) {
371 [ + - ]: 165882 : std::vector< std::vector< tk::real > > a( n.size() );
372 : 165882 : std::size_t j = 0;
373 [ + + ][ + - ]: 1769170 : for (auto i : n) a[ j++ ] = m_a[ tk::cref_find(m_lid,i) ];
[ + - ]
374 [ + - ][ + - ]: 165882 : thisProxy[ c ].comlim( std::vector<std::size_t>(begin(n),end(n)), a );
[ + - ]
375 : : }
376 : :
377 : 18283 : ownlim_complete();
378 : 18283 : }
379 : :
380 : : void
381 : 165882 : DistFCT::comlim( const std::vector< std::size_t >& gid,
382 : : const std::vector< std::vector< tk::real > >& A )
383 : : // *****************************************************************************
384 : : // Receive contributions of limited antidiffusive element contributions on
385 : : // chare-boundaries
386 : : //! \param[in] gid Global mesh node IDs at which we receive contributions
387 : : //! \param[in] A Partial contributions to antidiffusive element contributions to
388 : : //! chare-boundary nodes
389 : : //! \details This function receives contributions to m_a, which stores the
390 : : //! limited antidiffusive element contributions assembled to nodes (Lohner:
391 : : //! AEC^c), see also FluxCorrector::limit(). While m_a stores own
392 : : //! contributions, m_ac collects the neighbor chare contributions during
393 : : //! communication. This way work on m_a and m_ac is overlapped. The two are
394 : : //! combined in apply().
395 : : // *****************************************************************************
396 : : {
397 [ - + ][ - - ]: 165882 : Assert( A.size() == gid.size(), "Size mismatch" );
[ - - ][ - - ]
398 : :
399 : : using tk::operator+=;
400 : :
401 [ + + ]: 1769170 : for (std::size_t i=0; i<gid.size(); ++i) {
402 : 1603288 : auto bid = tk::cref_find( m_bid, gid[i] );
403 [ - + ][ - - ]: 1603288 : Assert( bid < m_ac.size(), "Indexing out of bounds" );
[ - - ][ - - ]
404 : 1603288 : m_ac[ bid ] += A[i];
405 : : }
406 : :
407 [ + + ]: 165882 : if (++m_nlim == m_nodeCommMap.size()) {
408 : 17784 : m_nlim = 0;
409 : 17784 : comlim_complete();
410 : : }
411 : 165882 : }
412 : :
413 : : void
414 : 18283 : DistFCT::apply()
415 : : // *****************************************************************************
416 : : // Apply limited antidiffusive element contributions
417 : : // *****************************************************************************
418 : : {
419 : : // Combine own and communicated contributions to A
420 [ + + ]: 1042826 : for (const auto& b : m_bid) {
421 [ + - ]: 1024543 : auto lid = tk::cref_find( m_lid, b.first );
422 : 1024543 : const auto& bac = m_ac[ b.second ];
423 [ + + ][ + - ]: 3131826 : for (ncomp_t c=0; c<m_a.nprop(); ++c) m_a(lid,c,0) += bac[c];
424 : : }
425 : :
426 : : // Update solution in host
427 [ + - ][ + - ]: 18283 : m_host[ thisIndex ].ckLocal()->update( m_a, std::move(m_dul) );
428 : 18283 : }
429 : :
430 : : std::tuple< std::vector< std::string >,
431 : : std::vector< std::vector< tk::real > >,
432 : : std::vector< std::string >,
433 : : std::vector< std::vector< tk::real > > >
434 : 2233 : DistFCT::fields() const
435 : : // *****************************************************************************
436 : : // Collect mesh output fields from FCT
437 : : //! \return Names and fields in mesh cells and nodes
438 : : // *****************************************************************************
439 : : {
440 [ + - ]: 2233 : auto fct_elemfields = m_fluxcorrector.fields( m_inpoel );
441 : :
442 : : using tuple_t = std::tuple< std::vector< std::string >,
443 : : std::vector< std::vector< tk::real > >,
444 : : std::vector< std::string >,
445 : : std::vector< std::vector< tk::real > > >;
446 : :
447 : 2233 : return tuple_t{ std::get< 0 >( fct_elemfields ),
448 : 2233 : std::get< 1 >( fct_elemfields ),
449 [ + - ]: 4466 : {}, {} };
450 : : }
451 : :
452 : : #include "NoWarning/distfct.def.h"
|
