Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/DiagCG.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 DiagCG for a PDE system with continuous Galerkin without a matrix
9 : : \details DiagCG advances a system of partial differential equations (PDEs)
10 : : using continuous Galerkin (CG) finite element (FE) spatial discretization
11 : : (using linear shapefunctions on tetrahedron elements) combined with a time
12 : : stepping scheme that is equivalent to the Lax-Wendroff (LW) scheme within
13 : : the unstructured-mesh FE context and treats discontinuities with
14 : : flux-corrected transport (FCT). Only the diagonal entries of the left-hand
15 : : side matrix are non-zero thus it does not need a matrix-based linear solver.
16 : : \see The documentation in DiagCG.h.
17 : : */
18 : : // *****************************************************************************
19 : :
20 : : #include "DiagCG.hpp"
21 : : #include "Vector.hpp"
22 : : #include "Reader.hpp"
23 : : #include "ContainerUtil.hpp"
24 : : #include "UnsMesh.hpp"
25 : : #include "Inciter/InputDeck/InputDeck.hpp"
26 : : #include "DerivedData.hpp"
27 : : #include "CGPDE.hpp"
28 : : #include "Discretization.hpp"
29 : : #include "DistFCT.hpp"
30 : : #include "DiagReducer.hpp"
31 : : #include "NodeBC.hpp"
32 : : #include "Refiner.hpp"
33 : : #include "Reorder.hpp"
34 : : #include "Integrate/Mass.hpp"
35 : : #include "FieldOutput.hpp"
36 : :
37 : : namespace inciter {
38 : :
39 : : extern ctr::InputDeck g_inputdeck;
40 : : extern ctr::InputDeck g_inputdeck_defaults;
41 : : extern std::vector< CGPDE > g_cgpde;
42 : :
43 : : } // inciter::
44 : :
45 : : using inciter::DiagCG;
46 : :
47 : 654 : DiagCG::DiagCG( const CProxy_Discretization& disc,
48 : : const std::map< int, std::vector< std::size_t > >& bface,
49 : : const std::map< int, std::vector< std::size_t > >& bnode,
50 : 654 : const std::vector< std::size_t >& triinpoel ) :
51 : : m_disc( disc ),
52 : : m_initial( 1 ),
53 : : m_nlhs( 0 ),
54 : : m_nrhs( 0 ),
55 : : m_nnorm( 0 ),
56 : : m_bnode( bnode ),
57 : : m_bface( bface ),
58 : : m_triinpoel( tk::remap(triinpoel,Disc()->Lid()) ),
59 : 1308 : m_u( Disc()->Gid().size(), g_inputdeck.get< tag::component >().nprop() ),
60 : : m_ul( m_u.nunk(), m_u.nprop() ),
61 : : m_du( m_u.nunk(), m_u.nprop() ),
62 [ + - ]: 654 : m_ue( Disc()->Inpoel().size()/4, m_u.nprop() ),
63 : : m_lhs( m_u.nunk(), m_u.nprop() ),
64 : : m_rhs( m_u.nunk(), m_u.nprop() ),
65 : : m_bcdir(),
66 : : m_lhsc(),
67 : : m_rhsc(),
68 : : m_difc(),
69 : : m_vol(),
70 : : m_bnorm(),
71 : : m_bnormc(),
72 : : m_symbcnodemap(),
73 : : m_symbcnodes(),
74 : : m_farfieldbcnodes(),
75 : : m_diag(),
76 : : m_boxnodes(),
77 [ + - ]: 654 : m_dtp( m_u.nunk(), 0.0 ),
78 : : m_tp( m_u.nunk(), g_inputdeck.get< tag::discr, tag::t0 >() ),
79 [ + - ][ + - ]: 1308 : m_finished( 0 )
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ]
80 : : // *****************************************************************************
81 : : // Constructor
82 : : //! \param[in] disc Discretization proxy
83 : : //! \param[in] bface Boundary-faces mapped to side set ids
84 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
85 : : //! \param[in] triinpoel Boundary-face connectivity
86 : : // *****************************************************************************
87 : : {
88 : 654 : usesAtSync = true; // enable migration at AtSync
89 : :
90 [ + - ]: 654 : auto d = Disc();
91 : :
92 : : // Perform optional operator-access-pattern mesh node reordering
93 [ + + ]: 654 : if (g_inputdeck.get< tag::discr, tag::operator_reorder >()) {
94 : :
95 : 47 : const auto& inpoel = d->Inpoel();
96 : :
97 : : // Create new local ids based on access pattern of PDE operators
98 : : std::unordered_map< std::size_t, std::size_t > map;
99 : : std::size_t n = 0;
100 : :
101 [ + + ]: 2237 : for (std::size_t e=0; e<inpoel.size()/4; ++e)
102 [ + + ]: 10950 : for (std::size_t i=0; i<4; ++i) {
103 : 8760 : std::size_t o = inpoel[e*4+i];
104 [ + + ][ + - ]: 8760 : if (map.find(o) == end(map)) map[o] = n++;
105 : : }
106 : :
107 : : Assert( map.size() == d->Gid().size(), "Map size mismatch" );
108 : :
109 : : // Remap data in bound Discretization object
110 [ + - ]: 47 : d->remap( map );
111 : : // Remap local ids in DistFCT
112 [ + - ][ + - ]: 47 : d->FCT()->remap( *d );
113 : : // Remap boundary triangle face connectivity
114 [ + - ]: 47 : tk::remap( m_triinpoel, map );
115 : :
116 : : }
117 : :
118 : : // Activate SDAG wait
119 [ + - ][ + - ]: 654 : thisProxy[ thisIndex ].wait4norm();
120 [ + - ][ + - ]: 654 : thisProxy[ thisIndex ].wait4lhs();
[ - - ]
121 : :
122 : : // Query nodes at which symmetry BCs are specified
123 [ + - ]: 654 : auto bn = d->bcnodes< tag::bc, tag::bcsym >( m_bface, m_triinpoel );
124 : : // Query nodes at which farfield BCs are specified
125 [ + - ]: 654 : auto far = d->bcnodes< tag::bc, tag::bcfarfield >( m_bface, m_triinpoel );
126 : :
127 : : // Merge BC data where boundary-point normals are required
128 [ - + ]: 654 : for (const auto& [s,n] : far) bn[s].insert( begin(n), end(n) );
129 : :
130 : : // Compute boundary point normals
131 [ + - ]: 654 : bnorm( bn );
132 : 654 : }
133 : :
134 : : void
135 : 654 : DiagCG::bnorm( const std::unordered_map< int,
136 : : std::unordered_set< std::size_t > >& bcnodes )
137 : : // *****************************************************************************
138 : : // Compute boundary point normals
139 : : //! \param[in] bcnodes Local node ids associated to side set ids at which BCs
140 : : //! are set that require normals
141 : : // *****************************************************************************
142 : : {
143 : 654 : auto d = Disc();
144 : :
145 [ + + ]: 1308 : m_bnorm = cg::bnorm( m_bface, m_triinpoel, d->Coord(), d->Gid(), bcnodes );
146 : :
147 : : // Send our nodal normal contributions to neighbor chares
148 [ + + ]: 654 : if (d->NodeCommMap().empty())
149 : 26 : comnorm_complete();
150 : : else
151 [ + + ]: 5952 : for (const auto& [ neighborchare, sharednodes ] : d->NodeCommMap()) {
152 : : std::unordered_map< int,
153 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > > exp;
154 [ + + ]: 111614 : for (auto i : sharednodes) {
155 [ + + ]: 110686 : for (const auto& [s,norms] : m_bnorm) {
156 : : auto j = norms.find(i);
157 [ + + ]: 5240 : if (j != end(norms)) exp[s][i] = j->second;
158 : : }
159 : : }
160 [ + - ][ + - ]: 10648 : thisProxy[ neighborchare ].comnorm( exp );
161 : : }
162 : :
163 : 654 : ownnorm_complete();
164 : 654 : }
165 : :
166 : : void
167 : 5324 : DiagCG::comnorm( const std::unordered_map< int,
168 : : std::unordered_map< std::size_t, std::array< tk::real, 4 > > >& innorm )
169 : : // *****************************************************************************
170 : : // Receive boundary point normals on chare-boundaries
171 : : //! \param[in] innorm Incoming partial sums of boundary point normal
172 : : //! contributions to normals (first 3 components), inverse distance squared
173 : : //! (4th component)
174 : : // *****************************************************************************
175 : : {
176 : : // Buffer up inccoming boundary-point normal vector contributions
177 [ + + ]: 5460 : for (const auto& [s,norms] : innorm) {
178 : : auto& bnorms = m_bnormc[s];
179 [ + + ]: 980 : for (const auto& [p,n] : norms) {
180 : : auto& bnorm = bnorms[p];
181 : 844 : bnorm[0] += n[0];
182 : 844 : bnorm[1] += n[1];
183 : 844 : bnorm[2] += n[2];
184 : 844 : bnorm[3] += n[3];
185 : : }
186 : : }
187 : :
188 [ + + ]: 5324 : if (++m_nnorm == Disc()->NodeCommMap().size()) {
189 : 628 : m_nnorm = 0;
190 : 628 : comnorm_complete();
191 : : }
192 : 5324 : }
193 : :
194 : : void
195 : 654 : DiagCG::normfinal()
196 : : // *****************************************************************************
197 : : // Finish computing boundary point normals
198 : : // *****************************************************************************
199 : : {
200 : 654 : auto d = Disc();
201 : 654 : const auto& lid = d->Lid();
202 : :
203 : : // Combine own and communicated contributions to boundary point normals
204 [ + + ]: 704 : for (const auto& [s,norms] : m_bnormc)
205 [ + + ]: 810 : for (const auto& [p,n] : norms) {
206 : : auto k = m_bnorm.find(s);
207 [ + + ]: 760 : if (k != end(m_bnorm)) {
208 : : auto j = k->second.find(p);
209 [ + - ]: 748 : if (j != end(k->second)) {
210 : : auto& norm = j->second;
211 : 748 : norm[0] += n[0];
212 : 748 : norm[1] += n[1];
213 : 748 : norm[2] += n[2];
214 : 748 : norm[3] += n[3];
215 : : }
216 : : }
217 : : }
218 : 654 : tk::destroy( m_bnormc );
219 : :
220 : : // Divie summed point normals by the sum of inverse distance squared
221 [ + + ]: 704 : for (auto& [s,norms] : m_bnorm)
222 [ + + ]: 2572 : for (auto& [p,n] : norms) {
223 : 2522 : n[0] /= n[3];
224 : 2522 : n[1] /= n[3];
225 : 2522 : n[2] /= n[3];
226 : : Assert( (n[0]*n[0] + n[1]*n[1] + n[2]*n[2] - 1.0) <
227 : : std::numeric_limits< tk::real >::epsilon(), "Non-unit normal" );
228 : : }
229 : :
230 : : // Replace global->local ids associated to boundary point normals
231 : : decltype(m_bnorm) bnorm;
232 [ + + ]: 704 : for (auto& [s,norms] : m_bnorm) {
233 : : auto& bnorms = bnorm[s];
234 [ + + ]: 2572 : for (auto&& [g,n] : norms)
235 : 2522 : bnorms[ tk::cref_find(lid,g) ] = std::move(n);
236 : : }
237 : : m_bnorm = std::move(bnorm);
238 : :
239 : : // Prepare unique set of symmetry BC nodes
240 [ + - ]: 1308 : m_symbcnodemap = d->bcnodes< tag::bc, tag::bcsym >( m_bface, m_triinpoel );
241 [ + + ]: 704 : for (const auto& [s,nodes] : m_symbcnodemap)
242 : : m_symbcnodes.insert( begin(nodes), end(nodes) );
243 : :
244 : : // Prepare unique set of farfield BC nodes
245 [ + - ]: 654 : auto far = d->bcnodes< tag::bc, tag::bcfarfield >( m_bface, m_triinpoel );
246 [ - + ]: 654 : for (const auto& [s,nodes] : far)
247 : : m_farfieldbcnodes.insert( begin(nodes), end(nodes) );
248 : :
249 : : // If farfield BC is set on a node, will not also set symmetry BC
250 [ - + ]: 654 : for (auto fn : m_farfieldbcnodes) {
251 : : m_symbcnodes.erase(fn);
252 [ - - ]: 0 : for (auto& [s,nodes] : m_symbcnodemap) nodes.erase(fn);
253 : : }
254 : :
255 : : // Signal the runtime system that the workers have been created
256 [ + - ]: 654 : std::vector< std::size_t > meshdata{ m_initial, d->MeshId() };
257 : 654 : contribute( meshdata, CkReduction::sum_ulong,
258 [ + - ][ + - ]: 1308 : CkCallback(CkReductionTarget(Transporter,comfinal), Disc()->Tr()) );
[ - - ]
259 : 654 : }
260 : :
261 : : void
262 : 666 : DiagCG::registerReducers()
263 : : // *****************************************************************************
264 : : // Configure Charm++ reduction types initiated from this chare array
265 : : //! \details Since this is a [initnode] routine, the runtime system executes the
266 : : //! routine exactly once on every logical node early on in the Charm++ init
267 : : //! sequence. Must be static as it is called without an object. See also:
268 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
269 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
270 : : // *****************************************************************************
271 : : {
272 : 666 : NodeDiagnostics::registerReducers();
273 : 666 : }
274 : :
275 : : void
276 : 9865 : DiagCG::ResumeFromSync()
277 : : // *****************************************************************************
278 : : // Return from migration
279 : : //! \details This is called when load balancing (LB) completes. The presence of
280 : : //! this function does not affect whether or not we block on LB.
281 : : // *****************************************************************************
282 : : {
283 [ - + ][ - - ]: 9865 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
284 : :
285 [ + - ]: 9865 : if (!g_inputdeck.get< tag::cmd, tag::nonblocking >()) next();
286 : 9865 : }
287 : :
288 : : void
289 : 654 : DiagCG::setup()
290 : : // *****************************************************************************
291 : : // Set and output initial conditions and mesh to file
292 : : // *****************************************************************************
293 : : {
294 : 654 : auto d = Disc();
295 : :
296 : : // Determine nodes inside user-defined IC box
297 [ + + ]: 1308 : for (auto& eq : g_cgpde) eq.IcBoxNodes( d->Coord(), m_boxnodes );
298 : :
299 : : // Compute volume of user-defined box IC
300 : 654 : d->boxvol( m_boxnodes );
301 : :
302 : : // Query time history field output labels from all PDEs integrated
303 : : const auto& hist_points = g_inputdeck.get< tag::history, tag::point >();
304 [ - + ]: 654 : if (!hist_points.empty()) {
305 : 0 : std::vector< std::string > histnames;
306 [ - - ]: 0 : for (const auto& eq : g_cgpde) {
307 : 0 : auto n = eq.histNames();
308 [ - - ]: 0 : histnames.insert( end(histnames), begin(n), end(n) );
309 : : }
310 [ - - ]: 0 : d->histheader( std::move(histnames) );
311 : : }
312 : 654 : }
313 : :
314 : : void
315 : 654 : DiagCG::box( tk::real v )
316 : : // *****************************************************************************
317 : : // Receive total box IC volume and set conditions in box
318 : : //! \param[in] v Total volume within user-specified box
319 : : // *****************************************************************************
320 : : {
321 : 654 : auto d = Disc();
322 : : const auto& coord = d->Coord();
323 : :
324 : : // Store user-defined box IC volume
325 : 654 : d->Boxvol() = v;
326 : :
327 : : // Set initial conditions for all PDEs
328 [ + + ]: 1308 : for (auto& eq : g_cgpde)
329 : 654 : eq.initialize( coord, m_u, d->T(), d->Boxvol(), m_boxnodes );
330 : :
331 : : // Apply symmetry BCs on initial conditions
332 [ + + ]: 1308 : for (const auto& eq : g_cgpde)
333 : 654 : eq.symbc( m_u, coord, m_bnorm, m_symbcnodes );
334 : : // Apply farfield BCs on initial conditions
335 [ + + ]: 1308 : for (const auto& eq : g_cgpde)
336 : 654 : eq.farfieldbc( m_u, coord, m_bnorm, m_farfieldbcnodes );
337 : :
338 : : // Output initial conditions to file (regardless of whether it was requested)
339 [ + - ][ + - ]: 1962 : writeFields( CkCallback(CkIndex_DiagCG::init(), thisProxy[thisIndex]) );
[ - + ][ - - ]
340 : 654 : }
341 : :
342 : : void
343 : 654 : DiagCG::init()
344 : : // *****************************************************************************
345 : : // Initially compute left hand side diagonal matrix
346 : : // *****************************************************************************
347 : : {
348 : 654 : lhs();
349 : 654 : }
350 : :
351 : : void
352 : 18283 : DiagCG::next()
353 : : // *****************************************************************************
354 : : // Continue to next time step
355 : : // *****************************************************************************
356 : : {
357 : 18283 : dt();
358 : 18283 : }
359 : :
360 : : void
361 : 730 : DiagCG::lhs()
362 : : // *****************************************************************************
363 : : // Compute the left-hand side of transport equations
364 : : // *****************************************************************************
365 : : {
366 : 730 : auto d = Disc();
367 : :
368 : : // Compute lumped mass lhs required for both high and low order solutions
369 [ + + ]: 1460 : m_lhs = tk::lump( m_u.nprop(), d->Coord(), d->Inpoel() );
370 : :
371 [ + + ]: 730 : if (d->NodeCommMap().empty())
372 : 28 : comlhs_complete();
373 : : else // send contributions of lhs to chare-boundary nodes to fellow chares
374 [ + + ]: 6498 : for (const auto& [c,n] : d->NodeCommMap()) {
375 [ + - ]: 5796 : std::vector< std::vector< tk::real > > l( n.size() );
376 : : std::size_t j = 0;
377 [ + + ][ + - ]: 258724 : for (auto i : n) l[ j++ ] = m_lhs[ tk::cref_find(d->Lid(),i) ];
[ - + ]
378 [ + - ][ + - ]: 11592 : thisProxy[c].comlhs( std::vector<std::size_t>(begin(n),end(n)), l );
[ + - ][ - + ]
379 : : }
380 : :
381 : 730 : ownlhs_complete();
382 : 730 : }
383 : :
384 : : void
385 : 5796 : DiagCG::comlhs( const std::vector< std::size_t >& gid,
386 : : const std::vector< std::vector< tk::real > >& L )
387 : : // *****************************************************************************
388 : : // Receive contributions to left-hand side diagonal matrix on chare-boundaries
389 : : //! \param[in] gid Global mesh node IDs at which we receive LHS contributions
390 : : //! \param[in] L Partial contributions of LHS to chare-boundary nodes
391 : : //! \details This function receives contributions to m_lhs, which stores the
392 : : //! diagonal (lumped) mass matrix at mesh nodes. While m_lhs stores
393 : : //! own contributions, m_lhsc collects the neighbor chare contributions during
394 : : //! communication. This way work on m_lhs and m_lhsc is overlapped. The two
395 : : //! are combined in lhsmerge().
396 : : // *****************************************************************************
397 : : {
398 : : Assert( L.size() == gid.size(), "Size mismatch" );
399 : :
400 : : using tk::operator+=;
401 : :
402 [ + + ]: 132260 : for (std::size_t i=0; i<gid.size(); ++i)
403 : 126464 : m_lhsc[ gid[i] ] += L[i];
404 : :
405 [ + + ]: 5796 : if (++m_nlhs == Disc()->NodeCommMap().size()) {
406 : 702 : m_nlhs = 0;
407 : 702 : comlhs_complete();
408 : : }
409 : 5796 : }
410 : :
411 : : void
412 : 730 : DiagCG::lhsmerge()
413 : : // *****************************************************************************
414 : : // The own and communication portion of the left-hand side is complete
415 : : // *****************************************************************************
416 : : {
417 : : // Combine own and communicated contributions to left hand side
418 [ + + ]: 91650 : for (const auto& b : m_lhsc) {
419 : 181840 : auto lid = tk::cref_find( Disc()->Lid(), b.first );
420 [ + + ]: 238716 : for (ncomp_t c=0; c<m_lhs.nprop(); ++c)
421 : 147796 : m_lhs(lid,c,0) += b.second[c];
422 : : }
423 : :
424 : : // Clear receive buffer
425 : 730 : tk::destroy(m_lhsc);
426 : :
427 : : // Continue after lhs is complete
428 [ + + ]: 730 : if (m_initial) {
429 : : // Start timer measuring time stepping wall clock time
430 : 654 : Disc()->Timer().zero();
431 : : // Zero grind-timer
432 : 654 : Disc()->grindZero();
433 : : // Continue to next time step
434 : 654 : next();
435 : : } else {
436 : 76 : lhs_complete();
437 : : }
438 : 730 : }
439 : :
440 : : void
441 : 18283 : DiagCG::dt()
442 : : // *****************************************************************************
443 : : // Compute time step size
444 : : // *****************************************************************************
445 : : {
446 : 18283 : tk::real mindt = std::numeric_limits< tk::real >::max();
447 : :
448 : 18283 : auto const_dt = g_inputdeck.get< tag::discr, tag::dt >();
449 : 18283 : auto def_const_dt = g_inputdeck_defaults.get< tag::discr, tag::dt >();
450 : : auto eps = std::numeric_limits< tk::real >::epsilon();
451 : :
452 : 18283 : auto d = Disc();
453 : :
454 : : // use constant dt if configured
455 [ + + ]: 18283 : if (std::abs(const_dt - def_const_dt) > eps) {
456 : :
457 : 3120 : mindt = const_dt;
458 : :
459 : : } else { // compute dt based on CFL
460 : :
461 : : // find the minimum dt across all PDEs integrated
462 [ + + ]: 30326 : for (const auto& eq : g_cgpde) {
463 : 15163 : auto eqdt = eq.dt( d->Coord(), d->Inpoel(), d->T(), d->Dtn(), m_u,
464 : : d->Vol(), d->Vol() );
465 [ + - ]: 15163 : if (eqdt < mindt) mindt = eqdt;
466 : : }
467 : :
468 : : }
469 : :
470 : : // Actiavate SDAG waits for time step
471 [ + - ]: 18283 : thisProxy[ thisIndex ].wait4rhs();
472 [ + - ]: 18283 : thisProxy[ thisIndex ].wait4out();
473 : :
474 : : // Activate SDAG-waits for FCT
475 : 18283 : d->FCT()->next();
476 : :
477 : : // Contribute to minimum dt across all chares the advance to next step
478 [ + - ]: 18283 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
479 : 18283 : CkCallback(CkReductionTarget(DiagCG,advance), thisProxy) );
480 : 18283 : }
481 : :
482 : : void
483 : 18283 : DiagCG::advance( tk::real newdt )
484 : : // *****************************************************************************
485 : : // Advance equations to next time step
486 : : //! \param[in] newdt Size of this new time step
487 : : // *****************************************************************************
488 : : {
489 : 18283 : auto d = Disc();
490 : :
491 : : // Set new time step size
492 : 18283 : d->setdt( newdt );
493 : :
494 : : // Compute rhs for next time step
495 : 18283 : rhs();
496 : 18283 : }
497 : :
498 : : void
499 : 18283 : DiagCG::rhs()
500 : : // *****************************************************************************
501 : : // Compute right-hand side of transport equations
502 : : // *****************************************************************************
503 : : {
504 : 18283 : auto d = Disc();
505 : 18283 : const auto& lid = d->Lid();
506 : 18283 : const auto& inpoel = d->Inpoel();
507 : :
508 : : // Sum nodal averages to elements (1st term of gather)
509 : : m_ue.fill( 0.0 );
510 [ + + ]: 9273135 : for (std::size_t e=0; e<inpoel.size()/4; ++e)
511 [ + + ]: 23143444 : for (ncomp_t c=0; c<m_u.nprop(); ++c)
512 [ + + ]: 69442960 : for (std::size_t a=0; a<4; ++a)
513 : 55554368 : m_ue(e,c,0) += m_u(inpoel[e*4+a],c,0)/4.0;
514 : :
515 : : // Scatter the right-hand side for chare-boundary cells only
516 : : m_rhs.fill( 0.0 );
517 [ + + ]: 36566 : for (const auto& eq : g_cgpde)
518 : 18283 : eq.rhs( d->T(), d->Dt(), d->Coord(), d->Inpoel(), m_u, m_ue, m_rhs );
519 : :
520 : : // Compute mass diffusion
521 : 18283 : auto dif = d->FCT()->diff( *d, m_u );
522 : :
523 : : // Query and match user-specified boundary conditions to side sets
524 [ + - ][ + + ]: 36566 : m_bcdir = match( m_u.nprop(), d->T(), d->Dt(), m_tp, m_dtp, d->Coord(),
525 [ + - ]: 18283 : lid, m_bnode, /* increment = */ true );
526 : :
527 : : // Send rhs data on chare-boundary nodes to fellow chares
528 [ + + ]: 18283 : if (d->NodeCommMap().empty())
529 [ + - ]: 499 : comrhs_complete();
530 : : else // send contributions of rhs to chare-boundary nodes to fellow chares
531 [ + + ]: 183666 : for (const auto& [c,n] : d->NodeCommMap()) {
532 [ + - ][ + - ]: 331764 : std::vector< std::vector< tk::real > > r( n.size() );
533 [ + - ]: 165882 : std::vector< std::vector< tk::real > > D( n.size() );
534 : : std::size_t j = 0;
535 [ + + ]: 1769170 : for (auto i : n) {
536 : 1603288 : auto k = tk::cref_find( lid, i );
537 [ + - ][ - + ]: 3206576 : r[j] = m_rhs[k];
[ + - ]
538 [ - + ]: 1603288 : D[j] = dif[k];
539 : 1603288 : ++j;
540 : : }
541 [ + - ][ + - ]: 331764 : thisProxy[c].comrhs( std::vector<std::size_t>(begin(n),end(n)), r, D );
[ + - ][ - + ]
542 : : }
543 : :
544 [ + - ]: 18283 : ownrhs_complete( dif );
545 : 18283 : }
546 : :
547 : : void
548 : 165882 : DiagCG::comrhs( const std::vector< std::size_t >& gid,
549 : : const std::vector< std::vector< tk::real > >& R,
550 : : const std::vector< std::vector< tk::real > >& D )
551 : : // *****************************************************************************
552 : : // Receive contributions to right-hand side vector on chare-boundaries
553 : : //! \param[in] gid Global mesh node IDs at which we receive RHS contributions
554 : : //! \param[in] R Partial contributions of RHS to chare-boundary nodes
555 : : //! \param[in] D Partial contributions to chare-boundary nodes
556 : : //! \details This function receives contributions to m_rhs, which stores the
557 : : //! right hand side vector at mesh nodes. While m_rhs stores own
558 : : //! contributions, m_rhsc collects the neighbor chare contributions during
559 : : //! communication. This way work on m_rhs and m_rhsc is overlapped. The two
560 : : //! are combined in solve(). This function also receives contributions to
561 : : //! mass diffusion term of the right hand side vector at mesh nodes.
562 : : // *****************************************************************************
563 : : {
564 : : Assert( R.size() == gid.size(), "Size mismatch" );
565 : : Assert( D.size() == gid.size(), "Size mismatch" );
566 : :
567 : : using tk::operator+=;
568 : :
569 [ + + ]: 1769170 : for (std::size_t i=0; i<gid.size(); ++i) {
570 : 1603288 : m_rhsc[ gid[i] ] += R[i];
571 : 1603288 : m_difc[ gid[i] ] += D[i];
572 : : }
573 : :
574 [ + + ]: 165882 : if (++m_nrhs == Disc()->NodeCommMap().size()) {
575 : 17784 : m_nrhs = 0;
576 : 17784 : comrhs_complete();
577 : : }
578 : 165882 : }
579 : :
580 : : void
581 : 18283 : DiagCG::solve( tk::Fields& dif )
582 : : // *****************************************************************************
583 : : // Solve low and high order diagonal systems
584 : : //! \param[in,out] dif Mass diffusion own contribution
585 : : // *****************************************************************************
586 : : {
587 : 18283 : const auto ncomp = m_rhs.nprop();
588 : :
589 : 18283 : auto d = Disc();
590 : :
591 : : // Combine own and communicated contributions to rhs
592 [ + + ]: 1042826 : for (const auto& b : m_rhsc) {
593 : 2049086 : auto lid = tk::cref_find( d->Lid(), b.first );
594 [ + + ]: 3131826 : for (ncomp_t c=0; c<ncomp; ++c) m_rhs(lid,c,0) += b.second[c];
595 : : }
596 : :
597 : : // Combine own and communicated contributions to mass diffusion
598 [ + + ]: 1042826 : for (const auto& b : m_difc) {
599 : 2049086 : auto lid = tk::cref_find( d->Lid(), b.first );
600 [ + + ]: 3131826 : for (ncomp_t c=0; c<ncomp; ++c) dif(lid,c,0) += b.second[c];
601 : : }
602 : :
603 : : // Clear receive buffers
604 : 18283 : tk::destroy(m_rhsc);
605 : 18283 : tk::destroy(m_difc);
606 : :
607 : : // Set Dirichlet BCs for lhs and both low and high order rhs vectors. Note
608 : : // that the low order rhs (more precisely the mass-diffusion term) is set to
609 : : // zero instead of the solution increment at Dirichlet BCs, because for the
610 : : // low order solution the right hand side is the sum of the high order right
611 : : // hand side and mass diffusion so the low order system is L = R + D, where L
612 : : // is the lumped mass matrix, R is the high order RHS, and D is
613 : : // mass diffusion, and R already will have the Dirichlet BC set.
614 [ + + ]: 528007 : for (const auto& [b,bc] : m_bcdir) {
615 [ + + ]: 2558044 : for (ncomp_t c=0; c<ncomp; ++c) {
616 [ + - ]: 2048320 : if (bc[c].first) {
617 : 2048320 : m_lhs( b, c, 0 ) = 1.0;
618 : 2048320 : m_rhs( b, c, 0 ) = bc[c].second;
619 : 2048320 : dif( b, c, 0 ) = 0.0;
620 : : }
621 : : }
622 : : }
623 : :
624 : : // Solve low and high order diagonal systems and update low order solution
625 [ + - ]: 18283 : auto dul = (m_rhs + dif) / m_lhs;
626 : :
627 [ + - ]: 18283 : m_ul = m_u + dul;
628 [ + - ]: 36566 : m_du = m_rhs / m_lhs;
629 : :
630 : : const auto& coord = d->Coord();
631 [ + + ]: 36566 : for (const auto& eq : g_cgpde) {
632 : : // Apply symmetry BCs
633 [ + - ]: 18283 : eq.symbc( dul, coord, m_bnorm, m_symbcnodes );
634 : : eq.symbc( m_ul, coord, m_bnorm, m_symbcnodes );
635 : : eq.symbc( m_du, coord, m_bnorm, m_symbcnodes );
636 : : // Apply farfield BCs
637 [ + - ]: 18283 : eq.farfieldbc( m_ul, coord, m_bnorm, m_farfieldbcnodes );
638 : : eq.farfieldbc( m_du, coord, m_bnorm, m_farfieldbcnodes );
639 : : }
640 : :
641 : : // Continue with FCT
642 [ + - ][ + - ]: 18283 : d->FCT()->aec( *d, m_du, m_u, m_bcdir, m_symbcnodemap, m_bnorm );
643 [ + - ][ + - ]: 18283 : d->FCT()->alw( m_u, m_ul, std::move(dul), thisProxy );
644 : 18283 : }
645 : :
646 : : void
647 : 2841 : DiagCG::writeFields( CkCallback c ) const
648 : : // *****************************************************************************
649 : : // Output mesh-based fields to file
650 : : //! \param[in] c Function to continue with after the write
651 : : // *****************************************************************************
652 : : {
653 [ + + ]: 2841 : if (g_inputdeck.get< tag::cmd, tag::benchmark >()) {
654 : :
655 : 608 : c.send();
656 : :
657 : : } else {
658 : :
659 : 2233 : auto d = Disc();
660 : : const auto& coord = d->Coord();
661 : :
662 : : // Query fields names requested by user
663 : 4466 : auto nodefieldnames = numericFieldNames( tk::Centering::NODE );
664 : : // Collect field output from numerical solution requested by user
665 [ + - ]: 6699 : auto nodefields = numericFieldOutput( m_u, tk::Centering::NODE );
666 : : // Collect field output names for analytical solutions
667 [ + + ]: 4466 : for (const auto& eq : g_cgpde)
668 [ + - ]: 2233 : analyticFieldNames( eq, tk::Centering::NODE, nodefieldnames );
669 : :
670 : : // Collect field output from analytical solutions (if exist)
671 : 2233 : auto t = d->T();
672 [ + + ]: 4466 : for (const auto& eq : g_cgpde)
673 [ + - ]: 2233 : analyticFieldOutput( eq, tk::Centering::NODE, coord[0], coord[1],
674 : : coord[2], t, nodefields );
675 : :
676 : : // Query and collect block and surface field names from PDEs integrated
677 : 2233 : std::vector< std::string > nodesurfnames;
678 [ + + ]: 4466 : for (const auto& eq : g_cgpde) {
679 : 2233 : auto s = eq.surfNames();
680 [ + - ]: 2233 : nodesurfnames.insert( end(nodesurfnames), begin(s), end(s) );
681 : : }
682 : :
683 : : // Collect node field solution
684 : : auto u = m_u;
685 : 2233 : std::vector< std::vector< tk::real > > nodesurfs;
686 [ + + ]: 4466 : for (const auto& eq : g_cgpde) {
687 [ + - ][ + - ]: 4466 : auto s = eq.surfOutput( tk::bfacenodes(m_bface,m_triinpoel), u );
688 [ + - ]: 2233 : nodesurfs.insert( end(nodesurfs), begin(s), end(s) );
689 : : }
690 : :
691 : : // Query refinement data
692 : : //auto dtref = g_inputdeck.get< tag::amr, tag::dtref >();
693 : :
694 : : std::tuple< std::vector< std::string >,
695 : : std::vector< std::vector< tk::real > >,
696 : : std::vector< std::string >,
697 : : std::vector< std::vector< tk::real > > > r;
698 [ + - ][ + - ]: 4466 : /*if (dtref)*/ r = d->Ref()->refinementFields();
[ + - ]
699 : :
700 : : auto& refinement_elemfieldnames = std::get< 0 >( r );
701 : : auto& refinement_elemfields = std::get< 1 >( r );
702 : : auto& refinement_nodefieldnames = std::get< 2 >( r );
703 : : auto& refinement_nodefields = std::get< 3 >( r );
704 : :
705 : : nodefieldnames.insert( end(nodefieldnames),
706 [ + - ][ + - ]: 2233 : begin(refinement_nodefieldnames), end(refinement_nodefieldnames) );
707 : : nodefields.insert( end(nodefields),
708 [ + - ][ + - ]: 2233 : begin(refinement_nodefields), end(refinement_nodefields) );
709 : :
710 : 2233 : auto elemfieldnames = std::move(refinement_elemfieldnames);
711 : 2233 : auto elemfields = std::move(refinement_elemfields);
712 : :
713 : : // Collect FCT field data (for debugging)
714 [ + - ][ + - ]: 2233 : auto f = d->FCT()->fields();
715 : :
716 : : const auto& fct_elemfieldnames = std::get< 0 >( f );
717 : : const auto& fct_elemfields = std::get< 1 >( f );
718 : : const auto& fct_nodefieldnames = std::get< 2 >( f );
719 : : const auto& fct_nodefields = std::get< 3 >( f );
720 : :
721 : : nodefieldnames.insert( end(nodefieldnames),
722 [ + - ][ + - ]: 2233 : begin(fct_nodefieldnames), end(fct_nodefieldnames) );
723 : : nodefields.insert( end(nodefields),
724 [ + - ][ + - ]: 2233 : begin(fct_nodefields), end(fct_nodefields) );
725 : :
726 : : elemfieldnames.insert( end(elemfieldnames),
727 [ + - ][ + - ]: 2233 : begin(fct_elemfieldnames), end(fct_elemfieldnames) );
728 : : elemfields.insert( end(elemfields),
729 [ + - ]: 2233 : begin(fct_elemfields), end(fct_elemfields) );
730 : :
731 : : Assert( elemfieldnames.size() == elemfields.size(), "Size mismatch" );
732 : : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
733 : :
734 : : // Send mesh and fields data (solution dump) for output to file
735 [ + - ]: 4466 : d->write( d->Inpoel(), coord, m_bface, tk::remap( m_bnode,d->Lid() ),
736 [ + - ]: 2233 : m_triinpoel, elemfieldnames, nodefieldnames, nodesurfnames,
737 : : elemfields, nodefields, nodesurfs, c );
738 : :
739 : : }
740 : 2841 : }
741 : :
742 : : void
743 : 18283 : DiagCG::update( const tk::Fields& a, [[maybe_unused]] tk::Fields&& dul )
744 : : // *****************************************************************************
745 : : // Prepare for next step
746 : : //! \param[in] a Limited antidiffusive element contributions
747 : : //! \param[in] dul Low order solution increment
748 : : // *****************************************************************************
749 : : {
750 : 18283 : auto d = Disc();
751 : :
752 : : // Verify that the change in the solution at those nodes where Dirichlet
753 : : // boundary conditions are set is exactly the amount the BCs prescribe
754 : : Assert( correctBC(a,dul,m_bcdir), "Dirichlet boundary condition incorrect" );
755 : :
756 : : // Apply limited antidiffusive element contributions to low order solution
757 : : auto un = m_u;
758 [ + + ]: 18283 : if (g_inputdeck.get< tag::discr, tag::fct >())
759 [ + - ]: 36546 : m_u = m_ul + a;
760 : : else
761 [ + - ]: 20 : m_u = m_u + m_du;
762 : :
763 : : // Compute diagnostics, e.g., residuals
764 : 18283 : auto diag_computed = m_diag.compute( *d, m_u, un, m_bnorm,
765 [ + - ]: 18283 : m_symbcnodes, m_farfieldbcnodes );
766 : : // Increase number of iterations and physical time
767 [ + - ]: 18283 : d->next();
768 : : // Continue to mesh refinement (if configured)
769 [ + + ][ + - ]: 18979 : if (!diag_computed) refine( std::vector< tk::real >( m_u.nprop(), 0.0 ) );
[ + - ][ - - ]
770 : 18283 : }
771 : :
772 : : void
773 : 18283 : DiagCG::refine( [[maybe_unused]] const std::vector< tk::real >& l2res )
774 : : // *****************************************************************************
775 : : // Optionally refine/derefine mesh
776 : : //! \param[in] l2res L2-norms of the residual for each scalar component
777 : : //! computed across the whole problem
778 : : // *****************************************************************************
779 : : {
780 : 18283 : auto d = Disc();
781 : :
782 : 18283 : auto dtref = g_inputdeck.get< tag::amr, tag::dtref >();
783 : 18283 : auto dtfreq = g_inputdeck.get< tag::amr, tag::dtfreq >();
784 : :
785 : : // if t>0 refinement enabled and we hit the dtref frequency
786 [ + + ][ + + ]: 18283 : if (dtref && !(d->It() % dtfreq)) { // h-refine
787 : :
788 : : // Activate SDAG waits for re-computing the left-hand side
789 [ + - ]: 76 : thisProxy[ thisIndex ].wait4lhs();
790 : :
791 : 76 : d->startvol();
792 [ + - ][ - + ]: 152 : d->Ref()->dtref( {}, m_bnode, {} );
[ - - ]
793 : 76 : d->refined() = 1;
794 : :
795 : : } else { // do not h-refine
796 : :
797 : 18207 : d->refined() = 0;
798 : 18207 : lhs_complete();
799 : 18207 : resized();
800 : :
801 : : }
802 : 18283 : }
803 : :
804 : : void
805 : 76 : DiagCG::resizePostAMR(
806 : : const std::vector< std::size_t >& /*ginpoel*/,
807 : : const tk::UnsMesh::Chunk& chunk,
808 : : const tk::UnsMesh::Coords& coord,
809 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
810 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
811 : : const std::set< std::size_t >& removedNodes,
812 : : const tk::NodeCommMap& nodeCommMap,
813 : : const std::map< int, std::vector< std::size_t > >& /*bface*/,
814 : : const std::map< int, std::vector< std::size_t > >& bnode,
815 : : const std::vector< std::size_t >& /*triinpoel*/ )
816 : : // *****************************************************************************
817 : : // Receive new mesh from Refiner
818 : : //! \param[in] ginpoel Mesh connectivity with global node ids
819 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
820 : : //! \param[in] coord New mesh node coordinates
821 : : //! \param[in] addedNodes Newly added mesh nodes and their parents (local ids)
822 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
823 : : //! \param[in] removedNodes Newly removed mesh nodes (local ids)
824 : : //! \param[in] nodeCommMap New node communication map
825 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
826 : : // *****************************************************************************
827 : : {
828 : 76 : auto d = Disc();
829 : :
830 : : // Set flag that indicates that we are during time stepping
831 : 76 : m_initial = 0;
832 : :
833 : : // Zero field output iteration count between two mesh refinement steps
834 : 76 : d->Itf() = 0;
835 : :
836 : : // Increase number of iterations with mesh refinement
837 : 76 : ++d->Itr();
838 : :
839 : : // Resize mesh data structures
840 : 76 : d->resizePostAMR( chunk, coord, nodeCommMap );
841 : :
842 : : // Remove newly removed nodes from solution vectors
843 : 76 : m_u.rm(removedNodes);
844 : 76 : m_ul.rm(removedNodes);
845 : 76 : m_du.rm(removedNodes);
846 : 76 : m_lhs.rm(removedNodes);
847 : 76 : m_rhs.rm(removedNodes);
848 : :
849 : : // Resize auxiliary solution vectors
850 : 76 : auto nelem = d->Inpoel().size()/4;
851 : 76 : auto npoin = coord[0].size();
852 : 76 : auto nprop = m_u.nprop();
853 : : m_u.resize( npoin );
854 : : m_ul.resize( npoin );
855 : : m_du.resize( npoin );
856 : : m_ue.resize( nelem );
857 : : m_lhs.resize( npoin );
858 : : m_rhs.resize( npoin );
859 : :
860 : : // Update solution on new mesh
861 [ + + ]: 50917 : for (const auto& n : addedNodes)
862 [ + + ]: 262690 : for (std::size_t c=0; c<nprop; ++c)
863 : 211849 : m_u(n.first,c,0) = (m_u(n.second[0],c,0) + m_u(n.second[1],c,0))/2.0;
864 : :
865 : : // Update physical-boundary node lists
866 : : m_bnode = bnode;
867 : :
868 : : // Resize FCT data structures
869 : 76 : d->FCT()->resize( npoin, nodeCommMap, d->Bid(), d->Lid(), d->Inpoel() );
870 : :
871 : 76 : auto meshid = d->MeshId();
872 [ + - ]: 152 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
873 : 76 : CkCallback(CkReductionTarget(Transporter,resized), d->Tr()) );
874 : 76 : }
875 : :
876 : : void
877 : 18283 : DiagCG::resized()
878 : : // *****************************************************************************
879 : : // Resizing data sutrctures after mesh refinement has been completed
880 : : // *****************************************************************************
881 : : {
882 : 18283 : resize_complete();
883 : 18283 : }
884 : :
885 : : void
886 : 18283 : DiagCG::out()
887 : : // *****************************************************************************
888 : : // Output mesh field data
889 : : // *****************************************************************************
890 : : {
891 : 18283 : auto d = Disc();
892 : :
893 : : // Output time history
894 [ + - ][ + - ]: 18283 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ]
895 : 0 : std::vector< std::vector< tk::real > > hist;
896 [ - - ]: 0 : for (const auto& eq : g_cgpde) {
897 [ - - ]: 0 : auto h = eq.histOutput( d->Hist(), d->Inpoel(), m_u );
898 [ - - ]: 0 : hist.insert( end(hist), begin(h), end(h) );
899 : : }
900 [ - - ]: 0 : d->history( std::move(hist) );
901 : : }
902 : :
903 : : // Output field data
904 [ + + ][ + - ]: 18283 : if (d->fielditer() or d->fieldtime() or d->fieldrange() or d->finished())
[ + - ][ + + ]
905 [ + - ][ + - ]: 6561 : writeFields( CkCallback(CkIndex_DiagCG::step(), thisProxy[thisIndex]) );
[ - + ][ - - ]
906 : : else
907 : 16096 : step();
908 : 18283 : }
909 : :
910 : : void
911 : 17629 : DiagCG::evalLB( int nrestart )
912 : : // *****************************************************************************
913 : : // Evaluate whether to do load balancing
914 : : //! \param[in] nrestart Number of times restarted
915 : : // *****************************************************************************
916 : : {
917 : 17629 : auto d = Disc();
918 : :
919 : : // Detect if just returned from a checkpoint and if so, zero timers
920 : 17629 : d->restarted( nrestart );
921 : :
922 : 17629 : const auto lbfreq = g_inputdeck.get< tag::cmd, tag::lbfreq >();
923 : 17629 : const auto nonblocking = g_inputdeck.get< tag::cmd, tag::nonblocking >();
924 : :
925 : : // Load balancing if user frequency is reached or after the second time-step
926 [ + + ][ + + ]: 17629 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
927 : :
928 : 9865 : AtSync();
929 [ - + ]: 9865 : if (nonblocking) next();
930 : :
931 : : } else {
932 : :
933 : 7764 : next();
934 : :
935 : : }
936 : 17629 : }
937 : :
938 : : void
939 : 17629 : DiagCG::evalRestart()
940 : : // *****************************************************************************
941 : : // Evaluate whether to save checkpoint/restart
942 : : // *****************************************************************************
943 : : {
944 : 17629 : auto d = Disc();
945 : :
946 : 17629 : const auto rsfreq = g_inputdeck.get< tag::cmd, tag::rsfreq >();
947 : 17629 : const auto benchmark = g_inputdeck.get< tag::cmd, tag::benchmark >();
948 : :
949 [ + + ][ - + ]: 17629 : if ( !benchmark && d->It() % rsfreq == 0 ) {
950 : :
951 : 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
952 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
953 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
954 : :
955 : : } else {
956 : :
957 : 17629 : evalLB( /* nrestart = */ -1 );
958 : :
959 : : }
960 : 17629 : }
961 : :
962 : : void
963 : 18283 : DiagCG::step()
964 : : // *****************************************************************************
965 : : // Evaluate whether to continue with next time step
966 : : // *****************************************************************************
967 : : {
968 : 18283 : auto d = Disc();
969 : :
970 : : // Output one-liner status report to screen
971 : 18283 : d->status();
972 : :
973 : 18283 : const auto term = g_inputdeck.get< tag::discr, tag::term >();
974 : 18283 : const auto nstep = g_inputdeck.get< tag::discr, tag::nstep >();
975 : : const auto eps = std::numeric_limits< tk::real >::epsilon();
976 : :
977 : : // If neither max iterations nor max time reached, continue, otherwise finish
978 [ + + ][ + + ]: 18283 : if (std::fabs(d->T()-term) > eps && d->It() < nstep) {
979 : :
980 : 17629 : evalRestart();
981 : :
982 : : } else {
983 : :
984 : 654 : auto meshid = d->MeshId();
985 [ + - ]: 1308 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
986 : 654 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
987 : :
988 : : }
989 : 18283 : }
990 : :
991 : : #include "NoWarning/diagcg.def.h"
|
