Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/ZalCG.cpp
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 ZalCG: Taylor-Galerkin, FCT, edge-based continuous Galerkin
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "XystBuildConfig.hpp"
14 : : #include "ZalCG.hpp"
15 : : #include "Vector.hpp"
16 : : #include "Reader.hpp"
17 : : #include "ContainerUtil.hpp"
18 : : #include "UnsMesh.hpp"
19 : : #include "ExodusIIMeshWriter.hpp"
20 : : #include "InciterConfig.hpp"
21 : : #include "DerivedData.hpp"
22 : : #include "Discretization.hpp"
23 : : #include "DiagReducer.hpp"
24 : : #include "IntegralReducer.hpp"
25 : : #include "Integrals.hpp"
26 : : #include "Refiner.hpp"
27 : : #include "Reorder.hpp"
28 : : #include "Around.hpp"
29 : : #include "Zalesak.hpp"
30 : : #include "Problems.hpp"
31 : : #include "EOS.hpp"
32 : : #include "BC.hpp"
33 : :
34 : : namespace inciter {
35 : :
36 : : extern ctr::Config g_cfg;
37 : :
38 : : static CkReduction::reducerType IntegralsMerger;
39 : :
40 : : } // inciter::
41 : :
42 : : using inciter::g_cfg;
43 : : using inciter::ZalCG;
44 : :
45 : 435 : ZalCG::ZalCG( const CProxy_Discretization& disc,
46 : : const std::map< int, std::vector< std::size_t > >& bface,
47 : : const std::map< int, std::vector< std::size_t > >& bnode,
48 : 435 : const std::vector< std::size_t >& triinpoel ) :
49 [ + - ]: 435 : m_disc( disc ),
50 : 435 : m_nrhs( 0 ),
51 : 435 : m_nnorm( 0 ),
52 : 435 : m_naec( 0 ),
53 : 435 : m_nalw( 0 ),
54 : 435 : m_nlim( 0 ),
55 : 435 : m_ndea( 0 ),
56 : 435 : m_nact( 0 ),
57 : 435 : m_todeactivate( 0 ),
58 : 435 : m_toreactivate( 0 ),
59 : 435 : m_deactivated( 0 ),
60 [ + - ]: 435 : m_bnode( bnode ),
61 [ + - ]: 435 : m_bface( bface ),
62 [ + - ][ + - ]: 435 : m_triinpoel( tk::remap( triinpoel, Disc()->Lid() ) ),
63 [ + - ][ + - ]: 435 : m_u( Disc()->Gid().size(), g_cfg.get< tag::problem_ncomp >() ),
64 [ + - ]: 435 : m_p( m_u.nunk(), m_u.nprop()*2 ),
65 [ + - ]: 435 : m_q( m_u.nunk(), m_u.nprop()*2 ),
66 [ + - ]: 435 : m_a( m_u.nunk(), m_u.nprop() ),
67 [ + - ]: 435 : m_rhs( m_u.nunk(), m_u.nprop() ),
68 [ + - ][ + - ]: 435 : m_vol( Disc()->Vol() ),
69 [ + - ]: 435 : m_dtp( m_u.nunk(), 0.0 ),
70 [ + - ]: 435 : m_tp( m_u.nunk(), g_cfg.get< tag::t0 >() ),
71 : 435 : m_finished( 0 ),
72 : 435 : m_freezeflow( 1.0 ),
73 : 3045 : m_fctfreeze( 0 )
74 : : // *****************************************************************************
75 : : // Constructor
76 : : //! \param[in] disc Discretization proxy
77 : : //! \param[in] bface Boundary-faces mapped to side sets used in the input file
78 : : //! \param[in] bnode Boundary-node lists mapped to side sets used in input file
79 : : //! \param[in] triinpoel Boundary-face connectivity where BCs set (global ids)
80 : : // *****************************************************************************
81 : : {
82 : 435 : usesAtSync = true; // enable migration at AtSync
83 : :
84 [ + - ]: 435 : auto d = Disc();
85 : :
86 : : // Compute total box IC volume
87 [ + - ]: 435 : d->boxvol();
88 : :
89 : : // Activate SDAG wait for initially computing integrals
90 [ + - ][ + - ]: 435 : thisProxy[ thisIndex ].wait4int();
91 : 435 : }
92 : :
93 : : void
94 : 435 : ZalCG::setupBC()
95 : : // *****************************************************************************
96 : : // Prepare boundary condition data structures
97 : : // *****************************************************************************
98 : : {
99 : : // Query Dirichlet BC nodes associated to side sets
100 : 435 : std::unordered_map< int, std::unordered_set< std::size_t > > dir;
101 [ + + ]: 615 : for (const auto& s : g_cfg.get< tag::bc_dir >()) {
102 [ + - ]: 180 : auto k = m_bface.find(s[0]);
103 [ + + ]: 180 : if (k != end(m_bface)) {
104 [ + - ]: 118 : auto& n = dir[ k->first ];
105 [ + + ]: 8268 : for (auto f : k->second) {
106 : 8150 : const auto t = m_triinpoel.data() + f*3;
107 [ + - ]: 8150 : n.insert( t[0] );
108 [ + - ]: 8150 : n.insert( t[1] );
109 [ + - ]: 8150 : n.insert( t[2] );
110 : : }
111 : : }
112 : : }
113 : :
114 : : // Augment Dirichlet BC nodes with nodes not necessarily part of faces
115 [ + - ]: 435 : const auto& lid = Disc()->Lid();
116 [ + + ]: 615 : for (const auto& s : g_cfg.get< tag::bc_dir >()) {
117 [ + - ]: 180 : auto k = m_bnode.find(s[0]);
118 [ + + ]: 180 : if (k != end(m_bnode)) {
119 [ + - ]: 118 : auto& n = dir[ k->first ];
120 [ + + ]: 13327 : for (auto g : k->second) {
121 [ + - ][ + - ]: 13209 : n.insert( tk::cref_find(lid,g) );
122 : : }
123 : : }
124 : : }
125 : :
126 : : // Collect unique set of nodes + Dirichlet BC components mask
127 : 435 : auto ncomp = m_u.nprop();
128 : 435 : auto nmask = ncomp + 1;
129 : 435 : const auto& dbc = g_cfg.get< tag::bc_dir >();
130 : 435 : std::unordered_map< std::size_t, std::vector< int > > dirbcset;
131 [ + + ]: 615 : for (const auto& mask : dbc) {
132 [ - + ][ - - ]: 180 : ErrChk( mask.size() == nmask, "Incorrect Dirichlet BC mask ncomp" );
[ - - ][ - - ]
133 [ + - ]: 180 : auto n = dir.find( mask[0] );
134 [ + + ]: 180 : if (n != end(dir))
135 [ + + ]: 5969 : for (auto p : n->second) {
136 [ + - ]: 5851 : auto& m = dirbcset[p];
137 [ + + ][ + - ]: 5851 : if (m.empty()) m.resize( ncomp, 0 );
138 [ + + ]: 40201 : for (std::size_t c=0; c<ncomp; ++c)
139 [ + + ]: 34350 : if (!m[c]) m[c] = mask[c+1]; // overwrite mask if 0 -> 1
140 : : }
141 : : }
142 : : // Compile streamable list of nodes + Dirichlet BC components mask
143 : 435 : tk::destroy( m_dirbcmasks );
144 [ + + ]: 5606 : for (const auto& [p,mask] : dirbcset) {
145 [ + - ]: 5171 : m_dirbcmasks.push_back( p );
146 [ + - ]: 5171 : m_dirbcmasks.insert( end(m_dirbcmasks), begin(mask), end(mask) );
147 : : }
148 [ - + ][ - - ]: 435 : ErrChk( m_dirbcmasks.size() % nmask == 0, "Dirichlet BC masks incomplete" );
[ - - ][ - - ]
149 : :
150 : : // Query pressure BC nodes associated to side sets
151 : 435 : const auto& pbc = g_cfg.get< tag::bc_pre >();
152 : 435 : const auto& pbc_sets = pbc.get< tag::sidesets >();
153 : 435 : std::unordered_map< int, std::unordered_set< std::size_t > > pre;
154 [ - + ]: 435 : for (const auto& ss : pbc_sets) {
155 [ - - ]: 0 : for (const auto& s : ss) {
156 [ - - ]: 0 : auto k = m_bface.find(s);
157 [ - - ]: 0 : if (k != end(m_bface)) {
158 [ - - ]: 0 : auto& n = pre[ k->first ];
159 [ - - ]: 0 : for (auto f : k->second) {
160 : 0 : const auto t = m_triinpoel.data() + f*3;
161 [ - - ]: 0 : n.insert( t[0] );
162 [ - - ]: 0 : n.insert( t[1] );
163 [ - - ]: 0 : n.insert( t[2] );
164 : : }
165 : : }
166 : : }
167 : : }
168 : :
169 : : // Augment Pressure BC nodes with nodes not necessarily part of faces
170 [ - + ]: 435 : for (const auto& s : pbc_sets) {
171 [ - - ]: 0 : auto k = m_bnode.find(s[0]);
172 [ - - ]: 0 : if (k != end(m_bnode)) {
173 [ - - ]: 0 : auto& n = pre[ k->first ];
174 [ - - ]: 0 : for (auto g : k->second) {
175 [ - - ][ - - ]: 0 : n.insert( tk::cref_find(lid,g) );
176 : : }
177 : : }
178 : : }
179 : :
180 : : // Prepare density and pressure values for pressure BC nodes
181 [ - + ]: 435 : if (!pbc_sets.empty()) {
182 : 0 : const auto& pbc_r = pbc.get< tag::density >();
183 [ - - ][ - - ]: 0 : ErrChk( pbc_r.size() == pbc_sets.size(), "Pressure BC density unspecified" );
[ - - ][ - - ]
184 : 0 : const auto& pbc_p = pbc.get< tag::pressure >();
185 [ - - ][ - - ]: 0 : ErrChk( pbc_p.size() == pbc_sets.size(), "Pressure BC pressure unspecified" );
[ - - ][ - - ]
186 : 0 : tk::destroy( m_prebcnodes );
187 : 0 : tk::destroy( m_prebcvals );
188 [ - - ]: 0 : for (const auto& [s,n] : pre) {
189 [ - - ]: 0 : m_prebcnodes.insert( end(m_prebcnodes), begin(n), end(n) );
190 [ - - ]: 0 : for (std::size_t p=0; p<pbc_sets.size(); ++p) {
191 [ - - ]: 0 : for (auto u : pbc_sets[p]) {
192 [ - - ]: 0 : if (s == u) {
193 [ - - ]: 0 : for (std::size_t i=0; i<n.size(); ++i) {
194 [ - - ]: 0 : m_prebcvals.push_back( pbc_r[p] );
195 [ - - ]: 0 : m_prebcvals.push_back( pbc_p[p] );
196 : : }
197 : : }
198 : : }
199 : : }
200 : : }
201 [ - - ][ - - ]: 0 : ErrChk( m_prebcnodes.size()*2 == m_prebcvals.size(),
[ - - ][ - - ]
202 : : "Pressure BC data incomplete" );
203 : : }
204 : :
205 : : // Query symmetry BC nodes associated to side sets
206 : 435 : std::unordered_map< int, std::unordered_set< std::size_t > > sym;
207 [ + + ]: 972 : for (auto s : g_cfg.get< tag::bc_sym >()) {
208 [ + - ]: 537 : auto k = m_bface.find(s);
209 [ + + ]: 537 : if (k != end(m_bface)) {
210 [ + - ]: 296 : auto& n = sym[ k->first ];
211 [ + + ]: 35707 : for (auto f : k->second) {
212 : 35411 : const auto t = m_triinpoel.data() + f*3;
213 [ + - ]: 35411 : n.insert( t[0] );
214 [ + - ]: 35411 : n.insert( t[1] );
215 [ + - ]: 35411 : n.insert( t[2] );
216 : : }
217 : : }
218 : : }
219 : :
220 : : // Query farfield BC nodes associated to side sets
221 : 435 : std::unordered_map< int, std::unordered_set< std::size_t > > far;
222 [ + + ]: 459 : for (auto s : g_cfg.get< tag::bc_far, tag::sidesets >()) {
223 [ + - ]: 24 : auto k = m_bface.find(s);
224 [ + + ]: 24 : if (k != end(m_bface)) {
225 [ + - ]: 7 : auto& n = far[ k->first ];
226 [ + + ]: 715 : for (auto f : k->second) {
227 : 708 : const auto t = m_triinpoel.data() + f*3;
228 [ + - ]: 708 : n.insert( t[0] );
229 [ + - ]: 708 : n.insert( t[1] );
230 [ + - ]: 708 : n.insert( t[2] );
231 : : }
232 : : }
233 : : }
234 : :
235 : : // Generate unique set of symmetry BC nodes
236 : 435 : tk::destroy( m_symbcnodeset );
237 [ + - ][ + + ]: 731 : for (const auto& [s,n] : sym) m_symbcnodeset.insert( begin(n), end(n) );
238 : : // Generate unique set of farfield BC nodes
239 : 435 : tk::destroy( m_farbcnodeset );
240 [ + - ][ + + ]: 442 : for (const auto& [s,n] : far) m_farbcnodeset.insert( begin(n), end(n) );
241 : :
242 : : // If farfield BC is set on a node, will not also set symmetry BC
243 [ + - ][ + + ]: 922 : for (auto i : m_farbcnodeset) m_symbcnodeset.erase(i);
244 : : // If pressure BC is set on a node, will not also set symmetry BC
245 [ - - ][ - + ]: 435 : for (auto i : m_prebcnodes) m_symbcnodeset.erase(i);
246 : 435 : }
247 : :
248 : : void
249 : 435 : ZalCG::feop()
250 : : // *****************************************************************************
251 : : // Start (re-)computing finite element domain and boundary operators
252 : : // *****************************************************************************
253 : : {
254 : 435 : auto d = Disc();
255 : :
256 : : // Prepare boundary conditions data structures
257 : 435 : setupBC();
258 : :
259 : : // Compute local contributions to boundary normals and integrals
260 : 435 : bndint();
261 : : // Compute contributions to domain edge integrals
262 : 435 : domint();
263 : :
264 : : // Send boundary point normal contributions to neighbor chares
265 [ + + ]: 435 : if (d->NodeCommMap().empty()) {
266 : 11 : comnorm_complete();
267 : : } else {
268 [ + + ]: 4574 : for (const auto& [c,nodes] : d->NodeCommMap()) {
269 : 4150 : decltype(m_bnorm) exp;
270 [ + + ]: 45118 : for (auto i : nodes) {
271 [ + + ]: 114322 : for (const auto& [s,b] : m_bnorm) {
272 [ + - ]: 73354 : auto k = b.find(i);
273 [ + + ][ + - ]: 73354 : if (k != end(b)) exp[s][i] = k->second;
[ + - ]
274 : : }
275 : : }
276 [ + - ][ + - ]: 4150 : thisProxy[c].comnorm( exp );
277 : 4150 : }
278 : : }
279 : 435 : ownnorm_complete();
280 : 435 : }
281 : :
282 : : void
283 : 435 : ZalCG::bndint()
284 : : // *****************************************************************************
285 : : //! Compute local contributions to boundary normals and integrals
286 : : // *****************************************************************************
287 : : {
288 [ + - ]: 435 : auto d = Disc();
289 : 435 : const auto& coord = d->Coord();
290 : 435 : const auto& gid = d->Gid();
291 : 435 : const auto& x = coord[0];
292 : 435 : const auto& y = coord[1];
293 : 435 : const auto& z = coord[2];
294 : :
295 : : // Lambda to compute the inverse distance squared between boundary face
296 : : // centroid and boundary point. Here p is the global node id and c is the
297 : : // the boundary face centroid.
298 : 171642 : auto invdistsq = [&]( const tk::real c[], std::size_t p ){
299 : 171642 : return 1.0 / ( (c[0] - x[p]) * (c[0] - x[p]) +
300 : 171642 : (c[1] - y[p]) * (c[1] - y[p]) +
301 : 171642 : (c[2] - z[p]) * (c[2] - z[p]) );
302 : 435 : };
303 : :
304 : 435 : tk::destroy( m_bnorm );
305 : 435 : tk::destroy( m_bndpoinint );
306 : 435 : tk::destroy( m_bndedgeint );
307 : :
308 [ + + ]: 1381 : for (const auto& [ setid, faceids ] : m_bface) { // for all side sets
309 [ + + ]: 58160 : for (auto f : faceids) { // for all side set triangles
310 : 57214 : const auto N = m_triinpoel.data() + f*3;
311 : : const std::array< tk::real, 3 >
312 : 57214 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
313 : 57214 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
314 : 57214 : auto n = tk::cross( ba, ca );
315 : 57214 : auto A2 = tk::length( n );
316 : 57214 : n[0] /= A2;
317 : 57214 : n[1] /= A2;
318 : 57214 : n[2] /= A2;
319 : : const tk::real centroid[3] = {
320 : 57214 : (x[N[0]] + x[N[1]] + x[N[2]]) / 3.0,
321 : 57214 : (y[N[0]] + y[N[1]] + y[N[2]]) / 3.0,
322 : 114428 : (z[N[0]] + z[N[1]] + z[N[2]]) / 3.0 };
323 : : // contribute all edges of triangle
324 [ + + ]: 228856 : for (const auto& [i,j] : tk::lpoet) {
325 : 171642 : auto p = N[i];
326 : 171642 : auto q = N[j];
327 : 171642 : tk::real r = invdistsq( centroid, p );
328 [ + - ]: 171642 : auto& v = m_bnorm[setid]; // associate side set id
329 [ + - ]: 171642 : auto& bn = v[gid[p]]; // associate global node id of bnd pnt
330 : 171642 : bn[0] += r * n[0]; // inv.dist.sq-weighted normal
331 : 171642 : bn[1] += r * n[1];
332 : 171642 : bn[2] += r * n[2];
333 : 171642 : bn[3] += r; // inv.dist.sq of node from centroid
334 [ + - ]: 171642 : auto& b = m_bndpoinint[gid[p]];// assoc global id of bnd point
335 : 171642 : b[0] += n[0] * A2 / 6.0; // bnd-point integral
336 : 171642 : b[1] += n[1] * A2 / 6.0;
337 : 171642 : b[2] += n[2] * A2 / 6.0;
338 : 171642 : tk::UnsMesh::Edge ed{ gid[p], gid[q] };
339 : 171642 : tk::real sig = 1.0;
340 [ + + ]: 171642 : if (ed[0] > ed[1]) {
341 : 85821 : std::swap( ed[0], ed[1] );
342 : 85821 : sig = -1.0;
343 : : }
344 [ + - ]: 171642 : auto& e = m_bndedgeint[ ed ];
345 : 171642 : e[0] += sig * n[0] * A2 / 24.0; // bnd-edge integral
346 : 171642 : e[1] += sig * n[1] * A2 / 24.0;
347 : 171642 : e[2] += sig * n[2] * A2 / 24.0;
348 : : }
349 : : }
350 : : }
351 : 435 : }
352 : :
353 : : void
354 : 435 : ZalCG::domint()
355 : : // *****************************************************************************
356 : : //! Compute local contributions to domain edge integrals
357 : : // *****************************************************************************
358 : : {
359 : 435 : auto d = Disc();
360 : :
361 : 435 : const auto& gid = d->Gid();
362 : 435 : const auto& inpoel = d->Inpoel();
363 : :
364 : 435 : const auto& coord = d->Coord();
365 : 435 : const auto& x = coord[0];
366 : 435 : const auto& y = coord[1];
367 : 435 : const auto& z = coord[2];
368 : :
369 : 435 : tk::destroy( m_domedgeint );
370 : :
371 [ + + ]: 393126 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
372 : 392691 : const auto N = inpoel.data() + e*4;
373 : : const std::array< tk::real, 3 >
374 : 392691 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
375 : 392691 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
376 : 392691 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
377 : 392691 : const auto J = tk::triple( ba, ca, da ); // J = 6V
378 [ - + ][ - - ]: 392691 : Assert( J > 0, "Element Jacobian non-positive" );
[ - - ][ - - ]
379 : : std::array< std::array< tk::real, 3 >, 4 > grad;
380 : 392691 : grad[1] = tk::cross( ca, da );
381 : 392691 : grad[2] = tk::cross( da, ba );
382 : 392691 : grad[3] = tk::cross( ba, ca );
383 [ + + ]: 1570764 : for (std::size_t i=0; i<3; ++i)
384 : 1178073 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
385 : 392691 : auto J120 = J/120.0;
386 [ + + ]: 2748837 : for (const auto& [p,q] : tk::lpoed) {
387 : 2356146 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
388 : 2356146 : tk::real sig = 1.0;
389 [ + + ]: 2356146 : if (ed[0] > ed[1]) {
390 : 908643 : std::swap( ed[0], ed[1] );
391 : 908643 : sig = -1.0;
392 : : }
393 [ + - ]: 2356146 : auto& n = m_domedgeint[ ed ];
394 : 2356146 : n[0] += sig * (grad[p][0] - grad[q][0]) / 48.0;
395 : 2356146 : n[1] += sig * (grad[p][1] - grad[q][1]) / 48.0;
396 : 2356146 : n[2] += sig * (grad[p][2] - grad[q][2]) / 48.0;
397 : 2356146 : n[3] += J120;
398 : : }
399 : : }
400 : 435 : }
401 : :
402 : : void
403 : 4150 : ZalCG::comnorm( const decltype(m_bnorm)& inbnd )
404 : : // *****************************************************************************
405 : : // Receive contributions to boundary point normals on chare-boundaries
406 : : //! \param[in] inbnd Incoming partial sums of boundary point normals
407 : : // *****************************************************************************
408 : : {
409 : : // Buffer up incoming boundary point normal vector contributions
410 [ + + ]: 7477 : for (const auto& [s,b] : inbnd) {
411 [ + - ]: 3327 : auto& bndnorm = m_bnormc[s];
412 [ + + ]: 16842 : for (const auto& [p,n] : b) {
413 [ + - ]: 13515 : auto& norm = bndnorm[p];
414 : 13515 : norm[0] += n[0];
415 : 13515 : norm[1] += n[1];
416 : 13515 : norm[2] += n[2];
417 : 13515 : norm[3] += n[3];
418 : : }
419 : : }
420 : :
421 [ + + ]: 4150 : if (++m_nnorm == Disc()->NodeCommMap().size()) {
422 : 424 : m_nnorm = 0;
423 : 424 : comnorm_complete();
424 : : }
425 : 4150 : }
426 : :
427 : : void
428 : 804 : ZalCG::registerReducers()
429 : : // *****************************************************************************
430 : : // Configure Charm++ reduction types initiated from this chare array
431 : : //! \details Since this is a [initnode] routine, the runtime system executes the
432 : : //! routine exactly once on every logical node early on in the Charm++ init
433 : : //! sequence. Must be static as it is called without an object. See also:
434 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
435 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
436 : : // *****************************************************************************
437 : : {
438 : 804 : NodeDiagnostics::registerReducers();
439 : 804 : IntegralsMerger = CkReduction::addReducer( integrals::mergeIntegrals );
440 : 804 : }
441 : :
442 : : void
443 : : // cppcheck-suppress unusedFunction
444 : 3614 : ZalCG::ResumeFromSync()
445 : : // *****************************************************************************
446 : : // Return from migration
447 : : //! \details This is called when load balancing (LB) completes. The presence of
448 : : //! this function does not affect whether or not we block on LB.
449 : : // *****************************************************************************
450 : : {
451 [ - + ][ - - ]: 3614 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
452 : :
453 [ + - ]: 3614 : if (!g_cfg.get< tag::nonblocking >()) dt();
454 : 3614 : }
455 : :
456 : : void
457 : 435 : ZalCG::setup( tk::real v )
458 : : // *****************************************************************************
459 : : // Start setup for solution
460 : : //! \param[in] v Total volume within user-specified box
461 : : // *****************************************************************************
462 : : {
463 : 435 : auto d = Disc();
464 : :
465 : : // Store user-defined box IC volume
466 : 435 : Disc()->Boxvol() = v;
467 : :
468 : : // Set initial conditions
469 : 435 : problems::initialize( d->Coord(), m_u, d->T(), /*meshid=*/0, d->BoxNodes() );
470 : :
471 : : // Query time history field output labels from all PDEs integrated
472 [ + + ]: 435 : if (!g_cfg.get< tag::histout, tag::points >().empty()) {
473 : : std::vector< std::string > var
474 [ + - ][ + + ]: 154 : {"density", "xvelocity", "yvelocity", "zvelocity", "energy", "pressure"};
[ - - ]
475 : 22 : auto ncomp = m_u.nprop();
476 [ - + ]: 22 : for (std::size_t c=5; c<ncomp; ++c)
477 [ - - ][ - - ]: 0 : var.push_back( "c" + std::to_string(c-5) );
[ - - ]
478 [ + - ]: 22 : d->histheader( std::move(var) );
479 : 22 : }
480 : :
481 : : // Compute finite element operators
482 : 435 : feop();
483 [ + - ][ + - ]: 479 : }
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
484 : :
485 : : void
486 : 435 : ZalCG::start()
487 : : // *****************************************************************************
488 : : // Start time stepping
489 : : // *****************************************************************************
490 : : {
491 : : // Set flag that indicates that we are now during time stepping
492 : 435 : Disc()->Initial( 0 );
493 : : // Start timer measuring time stepping wall clock time
494 : 435 : Disc()->Timer().zero();
495 : : // Zero grind-timer
496 : 435 : Disc()->grindZero();
497 : : // Continue to first time step
498 : 435 : dt();
499 : 435 : }
500 : :
501 : : void
502 : 435 : ZalCG::bnorm()
503 : : // *****************************************************************************
504 : : // Combine own and communicated portions of the boundary point normals
505 : : // *****************************************************************************
506 : : {
507 [ + - ]: 435 : const auto& lid = Disc()->Lid();
508 : :
509 : : // Combine own and communicated contributions to boundary point normals
510 [ + + ]: 1367 : for (const auto& [s,b] : m_bnormc) {
511 [ + - ]: 932 : auto& bndnorm = m_bnorm[s];
512 [ + + ]: 11628 : for (const auto& [g,n] : b) {
513 [ + - ]: 10696 : auto& norm = bndnorm[g];
514 : 10696 : norm[0] += n[0];
515 : 10696 : norm[1] += n[1];
516 : 10696 : norm[2] += n[2];
517 : 10696 : norm[3] += n[3];
518 : : }
519 : : }
520 : 435 : tk::destroy( m_bnormc );
521 : :
522 : : // Divide summed point normals by the sum of the inverse distance squared
523 [ + + ]: 1435 : for (auto& [s,b] : m_bnorm) {
524 [ + + ]: 39189 : for (auto& [g,n] : b) {
525 : 38189 : n[0] /= n[3];
526 : 38189 : n[1] /= n[3];
527 : 38189 : n[2] /= n[3];
528 [ - + ][ - - ]: 38189 : Assert( (n[0]*n[0] + n[1]*n[1] + n[2]*n[2] - 1.0) <
[ - - ][ - - ]
529 : : 1.0e+3*std::numeric_limits< tk::real >::epsilon(),
530 : : "Non-unit normal" );
531 : : }
532 : : }
533 : :
534 : : // Replace global->local ids associated to boundary point normals
535 : 435 : decltype(m_bnorm) loc;
536 [ + + ]: 1435 : for (auto& [s,b] : m_bnorm) {
537 [ + - ]: 1000 : auto& bnd = loc[s];
538 [ + + ]: 39189 : for (auto&& [g,n] : b) {
539 [ + - ][ + - ]: 38189 : bnd[ tk::cref_find(lid,g) ] = std::move(n);
540 : : }
541 : : }
542 : 435 : m_bnorm = std::move(loc);
543 : 435 : }
544 : :
545 : : void
546 : 435 : ZalCG::streamable()
547 : : // *****************************************************************************
548 : : // Convert integrals into streamable data structures
549 : : // *****************************************************************************
550 : : {
551 [ + - ]: 435 : const auto& lid = Disc()->Lid();
552 : :
553 : : // Convert boundary point integrals into streamable data structures
554 [ + - ]: 435 : m_bpoin.resize( m_bndpoinint.size() );
555 [ + - ]: 435 : m_bpint.resize( m_bndpoinint.size() * 3 );
556 : 435 : std::size_t i = 0;
557 [ + + ]: 34270 : for (const auto& [g,b] : m_bndpoinint) {
558 [ + - ]: 33835 : m_bpoin[i] = tk::cref_find( lid, g );
559 : 33835 : m_bpint[i*3+0] = b[0];
560 : 33835 : m_bpint[i*3+1] = b[1];
561 : 33835 : m_bpint[i*3+2] = b[2];
562 : 33835 : ++i;
563 : : }
564 : :
565 [ + - ]: 435 : m_besym.resize( m_triinpoel.size() );
566 : 435 : i = 0;
567 [ + + ]: 172077 : for (auto p : m_triinpoel) {
568 [ + - ]: 171642 : m_besym[i++] = static_cast< std::uint8_t >(m_symbcnodeset.count(p));
569 : : }
570 : :
571 : : // Query surface integral output nodes
572 : 435 : std::unordered_map< int, std::vector< std::size_t > > surfintnodes;
573 : 435 : const auto& is = g_cfg.get< tag::integout, tag::sidesets >();
574 [ + - ]: 435 : std::set< int > outsets( begin(is), end(is) );
575 [ - + ]: 435 : for (auto s : outsets) {
576 [ - - ]: 0 : auto m = m_bface.find(s);
577 [ - - ]: 0 : if (m != end(m_bface)) {
578 [ - - ]: 0 : auto& n = surfintnodes[ m->first ]; // associate set id
579 [ - - ]: 0 : for (auto f : m->second) { // face ids on side set
580 [ - - ]: 0 : n.push_back( m_triinpoel[f*3+0] ); // nodes on side set
581 [ - - ]: 0 : n.push_back( m_triinpoel[f*3+1] );
582 [ - - ]: 0 : n.push_back( m_triinpoel[f*3+2] );
583 : : }
584 : : }
585 : : }
586 [ - - ][ - + ]: 435 : for (auto& [s,n] : surfintnodes) tk::unique( n );
587 : : // Prepare surface integral data
588 : 435 : tk::destroy( m_surfint );
589 [ + - ]: 435 : const auto& gid = Disc()->Gid();
590 [ - + ]: 435 : for (auto&& [s,n] : surfintnodes) {
591 [ - - ]: 0 : auto& sint = m_surfint[s]; // associate set id
592 : 0 : auto& nodes = sint.first;
593 : 0 : auto& ndA = sint.second;
594 : 0 : nodes = std::move(n);
595 [ - - ]: 0 : ndA.resize( nodes.size()*3 );
596 : 0 : std::size_t a = 0;
597 [ - - ]: 0 : for (auto p : nodes) {
598 [ - - ]: 0 : const auto& b = tk::cref_find( m_bndpoinint, gid[p] );
599 : 0 : ndA[a*3+0] = b[0]; // store ni * dA
600 : 0 : ndA[a*3+1] = b[1];
601 : 0 : ndA[a*3+2] = b[2];
602 : 0 : ++a;
603 : : }
604 : : }
605 : 435 : tk::destroy( m_bndpoinint );
606 : :
607 : : // Generate edges along chare boundary
608 [ + - ]: 435 : chbnded();
609 : : // Adjust node volumes along inactive neighbor chares
610 [ + - ]: 435 : deavol();
611 : : // Generate superedges for domain integral
612 [ + - ]: 435 : domsuped();
613 : :
614 : : // Convert symmetry BC data to streamable data structures
615 : 435 : tk::destroy( m_symbcnodes );
616 : 435 : tk::destroy( m_symbcnorms );
617 [ + + ]: 20680 : for (auto p : m_symbcnodeset) {
618 [ + + ]: 61508 : for (const auto& s : g_cfg.get< tag::bc_sym >()) {
619 [ + - ]: 41263 : auto m = m_bnorm.find(s);
620 [ + + ]: 41263 : if (m != end(m_bnorm)) {
621 [ + - ]: 36196 : auto r = m->second.find(p);
622 [ + + ]: 36196 : if (r != end(m->second)) {
623 [ + - ]: 21521 : m_symbcnodes.push_back( p );
624 [ + - ]: 21521 : m_symbcnorms.push_back( r->second[0] );
625 [ + - ]: 21521 : m_symbcnorms.push_back( r->second[1] );
626 [ + - ]: 21521 : m_symbcnorms.push_back( r->second[2] );
627 : : }
628 : : }
629 : : }
630 : : }
631 : 435 : tk::destroy( m_symbcnodeset );
632 : :
633 : : // Convert farfield BC data to streamable data structures
634 : 435 : tk::destroy( m_farbcnodes );
635 : 435 : tk::destroy( m_farbcnorms );
636 [ + + ]: 922 : for (auto p : m_farbcnodeset) {
637 [ + + ]: 974 : for (const auto& s : g_cfg.get< tag::bc_far, tag::sidesets >()) {
638 [ + - ]: 487 : auto n = m_bnorm.find(s);
639 [ + - ]: 487 : if (n != end(m_bnorm)) {
640 [ + - ]: 487 : auto a = n->second.find(p);
641 [ + - ]: 487 : if (a != end(n->second)) {
642 [ + - ]: 487 : m_farbcnodes.push_back( p );
643 [ + - ]: 487 : m_farbcnorms.push_back( a->second[0] );
644 [ + - ]: 487 : m_farbcnorms.push_back( a->second[1] );
645 [ + - ]: 487 : m_farbcnorms.push_back( a->second[2] );
646 : : }
647 : : }
648 : : }
649 : : }
650 : 435 : tk::destroy( m_farbcnodeset );
651 : 435 : tk::destroy( m_bnorm );
652 : 435 : }
653 : :
654 : : void
655 : 435 : ZalCG::domsuped()
656 : : // *****************************************************************************
657 : : // Generate superedge-groups for domain-edge loops
658 : : //! \see See Lohner, Sec. 15.1.6.2, An Introduction to Applied CFD Techniques,
659 : : //! Wiley, 2008.
660 : : // *****************************************************************************
661 : : {
662 [ - + ][ - - ]: 435 : Assert( !m_domedgeint.empty(), "No domain edges to group" );
[ - - ][ - - ]
663 : :
664 : : #ifndef NDEBUG
665 : 435 : auto nedge = m_domedgeint.size();
666 : : #endif
667 : :
668 [ + - ]: 435 : const auto& inpoel = Disc()->Inpoel();
669 [ + - ]: 435 : const auto& lid = Disc()->Lid();
670 [ + - ]: 435 : const auto& gid = Disc()->Gid();
671 : :
672 : 435 : tk::destroy( m_dsupedge[0] );
673 : 435 : tk::destroy( m_dsupedge[1] );
674 : 435 : tk::destroy( m_dsupedge[2] );
675 : :
676 : 435 : tk::destroy( m_dsupint[0] );
677 : 435 : tk::destroy( m_dsupint[1] );
678 : 435 : tk::destroy( m_dsupint[2] );
679 : :
680 : 435 : tk::UnsMesh::FaceSet untri;
681 [ + + ]: 393126 : for (std::size_t e=0; e<inpoel.size()/4; e++) {
682 : 392691 : const auto N = inpoel.data() + e*4;
683 [ + - ][ + + ]: 1963455 : for (const auto& [a,b,c] : tk::lpofa) untri.insert( { N[a], N[b], N[c] } );
684 : : }
685 : :
686 [ + + ]: 393126 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
687 : 392691 : const auto N = inpoel.data() + e*4;
688 : 392691 : int f = 0;
689 : : tk::real sig[6];
690 [ + - ][ + + ]: 2748837 : decltype(m_domedgeint)::const_iterator d[6];
691 [ + + ]: 1158234 : for (const auto& [p,q] : tk::lpoed) {
692 : 1104399 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
693 [ + + ]: 1104399 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
694 [ + - ]: 1104399 : d[f] = m_domedgeint.find( ed );
695 [ + + ]: 1104399 : if (d[f] == end(m_domedgeint)) break; else ++f;
696 : : }
697 [ + + ]: 392691 : if (f == 6) {
698 [ + - ]: 53835 : m_dsupedge[0].push_back( N[0] );
699 [ + - ]: 53835 : m_dsupedge[0].push_back( N[1] );
700 [ + - ]: 53835 : m_dsupedge[0].push_back( N[2] );
701 [ + - ]: 53835 : m_dsupedge[0].push_back( N[3] );
702 [ + - ][ + + ]: 269175 : for (const auto& [a,b,c] : tk::lpofa) untri.erase( { N[a], N[b], N[c] } );
703 [ + + ]: 376845 : for (int ed=0; ed<6; ++ed) {
704 : 323010 : const auto& ded = d[ed]->second;
705 [ + - ]: 323010 : m_dsupint[0].push_back( sig[ed] * ded[0] );
706 [ + - ]: 323010 : m_dsupint[0].push_back( sig[ed] * ded[1] );
707 [ + - ]: 323010 : m_dsupint[0].push_back( sig[ed] * ded[2] );
708 [ + - ]: 323010 : m_dsupint[0].push_back( ded[3] );
709 [ + - ]: 323010 : m_domedgeint.erase( d[ed] );
710 : : }
711 : : }
712 : : }
713 : :
714 [ + + ]: 622713 : for (const auto& N : untri) {
715 : 622278 : int f = 0;
716 : : tk::real sig[3];
717 [ + - ][ + + ]: 2489112 : decltype(m_domedgeint)::const_iterator d[3];
718 [ + + ]: 985952 : for (const auto& [p,q] : tk::lpoet) {
719 : 955576 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
720 [ + + ]: 955576 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
721 [ + - ]: 955576 : d[f] = m_domedgeint.find( ed );
722 [ + + ]: 955576 : if (d[f] == end(m_domedgeint)) break; else ++f;
723 : : }
724 [ + + ]: 622278 : if (f == 3) {
725 [ + - ]: 30376 : m_dsupedge[1].push_back( N[0] );
726 [ + - ]: 30376 : m_dsupedge[1].push_back( N[1] );
727 [ + - ]: 30376 : m_dsupedge[1].push_back( N[2] );
728 [ + + ]: 121504 : for (int ed=0; ed<3; ++ed) {
729 : 91128 : const auto& ded = d[ed]->second;
730 [ + - ]: 91128 : m_dsupint[1].push_back( sig[ed] * ded[0] );
731 [ + - ]: 91128 : m_dsupint[1].push_back( sig[ed] * ded[1] );
732 [ + - ]: 91128 : m_dsupint[1].push_back( sig[ed] * ded[2] );
733 [ + - ]: 91128 : m_dsupint[1].push_back( ded[3] );
734 [ + - ]: 91128 : m_domedgeint.erase( d[ed] );
735 : : }
736 : : }
737 : : }
738 : :
739 [ + - ]: 435 : m_dsupedge[2].resize( m_domedgeint.size()*2 );
740 [ + - ]: 435 : m_dsupint[2].resize( m_domedgeint.size()*4 );
741 : 435 : std::size_t k = 0;
742 [ + + ]: 127344 : for (const auto& [ed,d] : m_domedgeint) {
743 : 126909 : auto e = m_dsupedge[2].data() + k*2;
744 [ + - ]: 126909 : e[0] = tk::cref_find( lid, ed[0] );
745 [ + - ]: 126909 : e[1] = tk::cref_find( lid, ed[1] );
746 : 126909 : auto n = m_dsupint[2].data() + k*4;
747 : 126909 : n[0] = d[0];
748 : 126909 : n[1] = d[1];
749 : 126909 : n[2] = d[2];
750 : 126909 : n[3] = d[3];
751 : 126909 : ++k;
752 : : }
753 : :
754 [ + - ]: 435 : if (g_cfg.get< tag::fct >()) {
755 : 435 : const auto ncomp = m_u.nprop();
756 [ + - ]: 435 : m_dsuplim[0].resize( m_dsupedge[0].size() * 6 * ncomp );
757 [ + - ]: 435 : m_dsuplim[1].resize( m_dsupedge[1].size() * 3 * ncomp );
758 [ + - ]: 435 : m_dsuplim[2].resize( m_dsupedge[2].size() * ncomp );
759 : : }
760 : :
761 : 435 : tk::destroy( m_domedgeint );
762 : :
763 : : //std::cout << std::setprecision(2)
764 : : // << "superedges: ntet:" << m_dsupedge[0].size()/4 << "(nedge:"
765 : : // << m_dsupedge[0].size()/4*6 << ","
766 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[0].size()/4*6 ) /
767 : : // static_cast< tk::real >( nedge )
768 : : // << "%) + ntri:" << m_dsupedge[1].size()/3
769 : : // << "(nedge:" << m_dsupedge[1].size() << ","
770 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[1].size() ) /
771 : : // static_cast< tk::real >( nedge )
772 : : // << "%) + nedge:"
773 : : // << m_dsupedge[2].size()/2 << "("
774 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[2].size()/2 ) /
775 : : // static_cast< tk::real >( nedge )
776 : : // << "%) = " << m_dsupedge[0].size()/4*6 + m_dsupedge[1].size() +
777 : : // m_dsupedge[2].size()/2 << " of "<< nedge << " total edges\n";
778 : :
779 [ - + ][ - - ]: 435 : Assert( m_dsupedge[0].size()/4*6 + m_dsupedge[1].size() +
[ - - ][ - - ]
780 : : m_dsupedge[2].size()/2 == nedge,
781 : : "Not all edges accounted for in superedge groups" );
782 : 435 : }
783 : :
784 : : void
785 : 435 : ZalCG::chbnded()
786 : : // *****************************************************************************
787 : : // Generate chare-boundary edge data structures for deactivation
788 : : // *****************************************************************************
789 : : {
790 [ + - ]: 435 : auto d = Disc();
791 [ + + ][ - + ]: 435 : if (not g_cfg.get< tag::deactivate >() or d->NodeCommMap().empty()) return;
[ + + ]
792 : :
793 [ + - ]: 19 : const auto& inpoel = Disc()->Inpoel();
794 [ + - ]: 19 : const auto& lid = Disc()->Lid();
795 [ + - ]: 19 : const auto& gid = Disc()->Gid();
796 : :
797 : : // Generate elements surrounding points
798 [ + - ]: 19 : auto esup = tk::genEsup( inpoel, 4 );
799 : :
800 : : // Collect edges of all tetrahedra surrounding chare-boundary points
801 : 19 : tk::destroy( m_chbndedge );
802 [ + + ]: 155 : for (const auto& [c,nodes] : d->NodeCommMap()) {
803 : : std::unordered_map< tk::UnsMesh::Edge, tk::real,
804 : 136 : tk::UnsMesh::Hash<2>, tk::UnsMesh::Eq<2> > edges;
805 [ + + ]: 8508 : for (auto g : nodes) {
806 [ + - ]: 8372 : auto i = tk::cref_find(lid,g);
807 [ + + ]: 89358 : for (auto e : tk::Around(esup,i)) {
808 : 80986 : const auto N = inpoel.data() + e*4;
809 [ + + ]: 566902 : for (const auto& [p,q] : tk::lpoed) {
810 : 485916 : tk::UnsMesh::Edge ged{ gid[N[p]], gid[N[q]] };
811 [ + - ][ + - ]: 485916 : edges[ { N[p], N[q] } ] = tk::cref_find(m_domedgeint,ged)[3];
812 : : }
813 : : }
814 : : }
815 [ + - ]: 136 : auto& ed = m_chbndedge[c];
816 [ + - ][ + + ]: 72203 : for (const auto& [e,sint] : edges) ed.emplace_back( e, sint );
817 : 136 : }
818 : 19 : }
819 : :
820 : : void
821 : 5765 : ZalCG::deavol()
822 : : // *****************************************************************************
823 : : // Adjust node volumes along inactive neighbor chares
824 : : // *****************************************************************************
825 : : {
826 : 5765 : auto d = Disc();
827 [ + + ][ - + ]: 5765 : if (not g_cfg.get< tag::deactivate >() or d->NodeCommMap().empty()) return;
[ + + ]
828 : :
829 : 259 : m_vol = d->Vol();
830 : 259 : const auto& cvolc = d->Cvolc();
831 [ + + ]: 944 : for (auto i : m_inactive) {
832 [ + - ][ + + ]: 33735 : for (const auto& [p,v] : tk::cref_find(cvolc,i)) {
833 : 33050 : m_vol[p] -= v;
834 : : }
835 : : }
836 : : }
837 : :
838 : : void
839 : : // cppcheck-suppress unusedFunction
840 : 435 : ZalCG::merge()
841 : : // *****************************************************************************
842 : : // Combine own and communicated portions of the integrals
843 : : // *****************************************************************************
844 : : {
845 : 435 : auto d = Disc();
846 : :
847 : : // Combine own and communicated contributions to boundary point normals
848 : 435 : bnorm();
849 : :
850 : : // Convert integrals into streamable data structures
851 : 435 : streamable();
852 : :
853 : : // Enforce boundary conditions using (re-)computed boundary data
854 : 435 : BC( m_u, d->T() );
855 : :
856 [ + - ]: 435 : if (d->Initial()) {
857 : : // Output initial conditions to file
858 [ + - ][ + - ]: 435 : writeFields( CkCallback(CkIndex_ZalCG::start(), thisProxy[thisIndex]) );
[ + - ][ + - ]
859 : : } else {
860 : 0 : feop_complete();
861 : : }
862 : 435 : }
863 : :
864 : : void
865 : 5765 : ZalCG::BC( tk::Fields& u, tk::real t )
866 : : // *****************************************************************************
867 : : // Apply boundary conditions
868 : : //! \param[in,out] u Solution to apply BCs to
869 : : //! \param[in] t Physical time
870 : : // *****************************************************************************
871 : : {
872 : 5765 : auto d = Disc();
873 : :
874 : : // Apply Dirichlet BCs
875 [ + - ]: 5765 : physics::dirbc( d->MeshId(), u, t, d->Coord(), d->BoxNodes(), m_dirbcmasks );
876 : :
877 : : // Apply symmetry BCs
878 : 5765 : physics::symbc( u, m_symbcnodes, m_symbcnorms, /*pos=*/1 );
879 : :
880 : : // Apply farfield BCs
881 : 5765 : physics::farbc( u, m_farbcnodes, m_farbcnorms );
882 : :
883 : : // Apply pressure BCs
884 : 5765 : physics::prebc( u, m_prebcnodes, m_prebcvals );
885 : 5765 : }
886 : :
887 : : void
888 : 6040 : ZalCG::dt()
889 : : // *****************************************************************************
890 : : // Compute time step size
891 : : // *****************************************************************************
892 : : {
893 : 6040 : tk::real mindt = std::numeric_limits< tk::real >::max();
894 : :
895 : 6040 : auto const_dt = g_cfg.get< tag::dt >();
896 : 6040 : auto eps = std::numeric_limits< tk::real >::epsilon();
897 [ + - ]: 6040 : auto d = Disc();
898 : :
899 [ + + ]: 6040 : if (not m_deactivated) {
900 : :
901 : : // Adjust node volumes along inactive neighbor chares for next step
902 [ + - ]: 5330 : deavol();
903 : :
904 : : // use constant dt if configured
905 [ + + ]: 5330 : if (std::abs(const_dt) > eps) {
906 : :
907 : : // cppcheck-suppress redundantInitialization
908 : 2920 : mindt = const_dt;
909 : :
910 : : } else {
911 : :
912 : 2410 : const auto& vol = d->Vol();
913 : 2410 : auto cfl = g_cfg.get< tag::cfl >();
914 : :
915 [ + + ]: 2410 : if (g_cfg.get< tag::steady >()) {
916 : :
917 [ + + ]: 792780 : for (std::size_t p=0; p<m_u.nunk(); ++p) {
918 [ + - ]: 792300 : auto r = m_u(p,0);
919 [ + - ]: 792300 : auto u = m_u(p,1)/r;
920 [ + - ]: 792300 : auto v = m_u(p,2)/r;
921 [ + - ]: 792300 : auto w = m_u(p,3)/r;
922 [ + - ]: 792300 : auto pr = eos::pressure( m_u(p,4) - 0.5*r*(u*u + v*v + w*w) );
923 : 792300 : auto c = eos::soundspeed( r, std::max(pr,0.0) );
924 : 792300 : auto L = std::cbrt( vol[p] );
925 : 792300 : auto vel = std::sqrt( u*u + v*v + w*w );
926 : 792300 : m_dtp[p] = L / std::max( vel+c, 1.0e-8 ) * cfl;
927 : : }
928 [ + - ]: 480 : mindt = *std::min_element( begin(m_dtp), end(m_dtp) );
929 : :
930 : : } else {
931 : :
932 [ + + ]: 832470 : for (std::size_t p=0; p<m_u.nunk(); ++p) {
933 [ + - ]: 830540 : auto r = m_u(p,0);
934 [ + - ]: 830540 : auto u = m_u(p,1)/r;
935 [ + - ]: 830540 : auto v = m_u(p,2)/r;
936 [ + - ]: 830540 : auto w = m_u(p,3)/r;
937 [ + - ]: 830540 : auto pr = eos::pressure( m_u(p,4) - 0.5*r*(u*u + v*v + w*w) );
938 : 830540 : auto c = eos::soundspeed( r, std::max(pr,0.0) );
939 : 830540 : auto L = std::cbrt( vol[p] );
940 : 830540 : auto vel = std::sqrt( u*u + v*v + w*w );
941 : 830540 : auto euler_dt = L / std::max( vel+c, 1.0e-8 );
942 : 830540 : mindt = std::min( mindt, euler_dt );
943 : : }
944 : 1930 : mindt *= cfl;
945 : :
946 : : }
947 : :
948 : : // Freeze flow if configured and apply multiplier on scalar(s)
949 [ + + ]: 2410 : if (d->T() > g_cfg.get< tag::freezetime >()) {
950 : 2267 : m_freezeflow = g_cfg.get< tag::freezeflow >();
951 : : }
952 : 2410 : mindt *= m_freezeflow;
953 : :
954 : : }
955 : :
956 : : }
957 : :
958 : : // Actiavate SDAG waits for next time step
959 [ + - ][ + - ]: 6040 : thisProxy[ thisIndex ].wait4rhs();
960 [ + - ][ + - ]: 6040 : thisProxy[ thisIndex ].wait4aec();
961 [ + - ][ + - ]: 6040 : thisProxy[ thisIndex ].wait4alw();
962 [ + - ][ + - ]: 6040 : thisProxy[ thisIndex ].wait4sol();
963 [ + - ][ + - ]: 6040 : thisProxy[ thisIndex ].wait4dea();
964 [ + - ][ + - ]: 6040 : thisProxy[ thisIndex ].wait4act();
965 [ + - ][ + - ]: 6040 : thisProxy[ thisIndex ].wait4step();
966 : :
967 : : // Contribute to minimum dt across all chares and advance to next step
968 [ + - ]: 6040 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
969 [ + - ][ + - ]: 12080 : CkCallback(CkReductionTarget(ZalCG,advance), thisProxy) );
970 : 6040 : }
971 : :
972 : : void
973 : 6040 : ZalCG::advance( tk::real newdt )
974 : : // *****************************************************************************
975 : : // Advance equations to next time step
976 : : //! \param[in] newdt The smallest dt across the whole problem
977 : : // *****************************************************************************
978 : : {
979 : : // Detect blowup
980 : 6040 : auto eps = std::numeric_limits< tk::real >::epsilon();
981 [ - + ]: 6040 : if (newdt < eps) m_finished = 1;
982 : :
983 : : // Set new time step size
984 : 6040 : Disc()->setdt( newdt );
985 : :
986 [ + + ]: 6040 : if (m_deactivated) solve(); else rhs();
987 : 6040 : }
988 : :
989 : : void
990 : 5330 : ZalCG::rhs()
991 : : // *****************************************************************************
992 : : // Compute right-hand side of transport equations
993 : : // *****************************************************************************
994 : : {
995 : 5330 : auto d = Disc();
996 : 5330 : const auto& lid = d->Lid();
997 : :
998 : : // Compute own portion of right-hand side for all equations
999 : 5330 : zalesak::rhs( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_besym,
1000 : 5330 : d->T(), d->Dt(), m_tp, m_dtp, m_u, m_rhs );
1001 : :
1002 : : // Communicate rhs to other chares on chare-boundary
1003 [ + + ][ + + ]: 5330 : if (d->NodeCommMap().empty() or m_inactive.size() == d->NodeCommMap().size()){
[ + + ]
1004 : 145 : comrhs_complete();
1005 : : } else {
1006 [ + + ]: 51830 : for (const auto& [c,n] : d->NodeCommMap()) {
1007 [ + - ][ + + ]: 46645 : if (m_inactive.count(c)) continue;
1008 : 45980 : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1009 [ + - ][ + - ]: 629110 : for (auto g : n) exp[g] = m_rhs[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1010 [ + - ][ + - ]: 45980 : thisProxy[c].comrhs( exp );
1011 : 45980 : }
1012 : : }
1013 : 5330 : ownrhs_complete();
1014 : 5330 : }
1015 : :
1016 : : void
1017 : 45980 : ZalCG::comrhs(
1018 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& inrhs )
1019 : : // *****************************************************************************
1020 : : // Receive contributions to right-hand side vector on chare-boundaries
1021 : : //! \param[in] inrhs Partial contributions of RHS to chare-boundary nodes. Key:
1022 : : //! global mesh node IDs, value: contributions for all scalar components.
1023 : : // *****************************************************************************
1024 : : {
1025 : : using tk::operator+=;
1026 [ + - ][ + - ]: 629110 : for (const auto& [g,r] : inrhs) m_rhsc[g] += r;
[ + + ]
1027 : :
1028 : : // When we have heard from all chares we communicate with, this chare is done
1029 [ + + ]: 45980 : if (++m_nrhs + m_inactive.size() == Disc()->NodeCommMap().size()) {
1030 : 5185 : m_nrhs = 0;
1031 : 5185 : comrhs_complete();
1032 : : }
1033 : 45980 : }
1034 : :
1035 : : void
1036 : 5330 : ZalCG::fct()
1037 : : // *****************************************************************************
1038 : : // Continue with flux-corrected transport if enabled
1039 : : // *****************************************************************************
1040 : : {
1041 : 5330 : auto d = Disc();
1042 : 5330 : const auto& lid = d->Lid();
1043 : :
1044 : : // Combine own and communicated contributions to rhs
1045 [ + + ]: 442280 : for (const auto& [g,r] : m_rhsc) {
1046 [ + - ]: 436950 : auto i = tk::cref_find( lid, g );
1047 [ + - ][ + + ]: 2656840 : for (std::size_t c=0; c<r.size(); ++c) m_rhs(i,c) += r[c];
1048 : : }
1049 : 5330 : tk::destroy(m_rhsc);
1050 : :
1051 [ + - ]: 5330 : if (g_cfg.get< tag::fct >()) aec(); else solve();
1052 : 5330 : }
1053 : :
1054 : : void
1055 : : // cppcheck-suppress unusedFunction
1056 : 5330 : ZalCG::aec()
1057 : : // *****************************************************************************
1058 : : // Compute antidiffusive contributions: P+/-
1059 : : // *****************************************************************************
1060 : : {
1061 : 5330 : auto d = Disc();
1062 : 5330 : const auto ncomp = m_u.nprop();
1063 : 5330 : const auto& lid = d->Lid();
1064 : :
1065 : : // Antidiffusive contributions: P+/-
1066 : :
1067 : 5330 : auto ctau = g_cfg.get< tag::fctdif >();
1068 : 5330 : m_p.fill( 0.0 );
1069 : :
1070 : : // tetrahedron superedges
1071 [ + + ]: 984945 : for (std::size_t e=0; e<m_dsupedge[0].size()/4; ++e) {
1072 : 979615 : const auto N = m_dsupedge[0].data() + e*4;
1073 : 979615 : const auto D = m_dsupint[0].data();
1074 : 979615 : std::size_t i = 0;
1075 [ + + ]: 6857305 : for (const auto& [p,q] : tk::lpoed) {
1076 : 5877690 : auto dif = D[(e*6+i)*4+3];
1077 [ + + ]: 35463900 : for (std::size_t c=0; c<ncomp; ++c) {
1078 [ + - ][ + - ]: 29586210 : auto aec = -dif * ctau * (m_u(N[p],c) - m_u(N[q],c));
1079 : 29586210 : auto a = c*2;
1080 : 29586210 : auto b = a+1;
1081 [ + + ]: 29586210 : if (aec > 0.0) std::swap(a,b);
1082 [ + - ]: 29586210 : m_p(N[p],a) -= aec;
1083 [ + - ]: 29586210 : m_p(N[q],b) += aec;
1084 : : }
1085 : 5877690 : ++i;
1086 : : }
1087 : : }
1088 : :
1089 : : // triangle superedges
1090 [ + + ]: 533300 : for (std::size_t e=0; e<m_dsupedge[1].size()/3; ++e) {
1091 : 527970 : const auto N = m_dsupedge[1].data() + e*3;
1092 : 527970 : const auto D = m_dsupint[1].data();
1093 : 527970 : std::size_t i = 0;
1094 [ + + ]: 2111880 : for (const auto& [p,q] : tk::lpoet) {
1095 : 1583910 : auto dif = D[(e*3+i)*4+3];
1096 [ + + ]: 9597360 : for (std::size_t c=0; c<ncomp; ++c) {
1097 [ + - ][ + - ]: 8013450 : auto aec = -dif * ctau * (m_u(N[p],c) - m_u(N[q],c));
1098 : 8013450 : auto a = c*2;
1099 : 8013450 : auto b = a+1;
1100 [ + + ]: 8013450 : if (aec > 0.0) std::swap(a,b);
1101 [ + - ]: 8013450 : m_p(N[p],a) -= aec;
1102 [ + - ]: 8013450 : m_p(N[q],b) += aec;
1103 : : }
1104 : 1583910 : ++i;
1105 : : }
1106 : : }
1107 : :
1108 : : // edges
1109 [ + + ]: 2272950 : for (std::size_t e=0; e<m_dsupedge[2].size()/2; ++e) {
1110 : 2267620 : const auto N = m_dsupedge[2].data() + e*2;
1111 : 2267620 : const auto dif = m_dsupint[2][e*4+3];
1112 [ + + ]: 13712040 : for (std::size_t c=0; c<ncomp; ++c) {
1113 [ + - ][ + - ]: 11444420 : auto aec = -dif * ctau * (m_u(N[0],c) - m_u(N[1],c));
1114 : 11444420 : auto a = c*2;
1115 : 11444420 : auto b = a+1;
1116 [ + + ]: 11444420 : if (aec > 0.0) std::swap(a,b);
1117 [ + - ]: 11444420 : m_p(N[0],a) -= aec;
1118 [ + - ]: 11444420 : m_p(N[1],b) += aec;
1119 : : }
1120 : : }
1121 : :
1122 : : // Apply symmetry BCs on AEC
1123 [ + + ]: 376690 : for (std::size_t i=0; i<m_symbcnodes.size(); ++i) {
1124 : 371360 : auto p = m_symbcnodes[i];
1125 : 371360 : auto nx = m_symbcnorms[i*3+0];
1126 : 371360 : auto ny = m_symbcnorms[i*3+1];
1127 : 371360 : auto nz = m_symbcnorms[i*3+2];
1128 : 371360 : auto rvnp = m_p(p,2)*nx + m_p(p,4)*ny + m_p(p,6)*nz;
1129 : 371360 : auto rvnn = m_p(p,3)*nx + m_p(p,5)*ny + m_p(p,7)*nz;
1130 : 371360 : m_p(p,2) -= rvnp * nx;
1131 : 371360 : m_p(p,3) -= rvnn * nx;
1132 : 371360 : m_p(p,4) -= rvnp * ny;
1133 : 371360 : m_p(p,5) -= rvnn * ny;
1134 : 371360 : m_p(p,6) -= rvnp * nz;
1135 : 371360 : m_p(p,7) -= rvnn * nz;
1136 : : }
1137 : :
1138 : : // Communicate antidiffusive edge and low-order solution contributions
1139 [ + + ][ + + ]: 5330 : if (d->NodeCommMap().empty() or m_inactive.size() == d->NodeCommMap().size()){
[ + + ]
1140 : 145 : comaec_complete();
1141 : : } else {
1142 [ + + ]: 51830 : for (const auto& [c,n] : d->NodeCommMap()) {
1143 [ + - ][ + + ]: 46645 : if (m_inactive.count(c)) continue;
1144 : 45980 : decltype(m_pc) exp;
1145 [ + - ][ + - ]: 629110 : for (auto g : n) exp[g] = m_p[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1146 [ + - ][ + - ]: 45980 : thisProxy[c].comaec( exp );
1147 : 45980 : }
1148 : : }
1149 : 5330 : ownaec_complete();
1150 : 5330 : }
1151 : :
1152 : : void
1153 : 45980 : ZalCG::comaec( const std::unordered_map< std::size_t,
1154 : : std::vector< tk::real > >& inaec )
1155 : : // *****************************************************************************
1156 : : // Receive antidiffusive and low-order contributions on chare-boundaries
1157 : : //! \param[in] inaec Partial contributions of antidiffusive edge and low-order
1158 : : //! solution contributions on chare-boundary nodes. Key: global mesh node IDs,
1159 : : //! value: 0: antidiffusive contributions, 1: low-order solution.
1160 : : // *****************************************************************************
1161 : : {
1162 : : using tk::operator+=;
1163 [ + - ][ + - ]: 629110 : for (const auto& [g,a] : inaec) m_pc[g] += a;
[ + + ]
1164 : :
1165 : : // When we have heard from all chares we communicate with, this chare is done
1166 [ + + ]: 45980 : if (++m_naec + m_inactive.size() == Disc()->NodeCommMap().size()) {
1167 : 5185 : m_naec = 0;
1168 : 5185 : comaec_complete();
1169 : : }
1170 : 45980 : }
1171 : :
1172 : : void
1173 : 5330 : ZalCG::alw()
1174 : : // *****************************************************************************
1175 : : // Compute allowed limits, Q+/-
1176 : : // *****************************************************************************
1177 : : {
1178 : 5330 : auto d = Disc();
1179 : 5330 : const auto steady = g_cfg.get< tag::steady >();
1180 : 5330 : const auto npoin = m_u.nunk();
1181 : 5330 : const auto ncomp = m_u.nprop();
1182 : 5330 : const auto& lid = d->Lid();
1183 : :
1184 : : // Combine own and communicated contributions to antidiffusive contributions
1185 : : // and low-order solution
1186 [ + + ]: 442280 : for (const auto& [g,p] : m_pc) {
1187 [ + - ]: 436950 : auto i = tk::cref_find( lid, g );
1188 [ + - ][ + + ]: 4876730 : for (std::size_t c=0; c<p.size(); ++c) m_p(i,c) += p[c];
1189 : : }
1190 : 5330 : tk::destroy(m_pc);
1191 : :
1192 : : // Finish computing antidiffusive contributions and low-order solution
1193 : 5330 : auto dt = d->Dt();
1194 [ + + ]: 1705150 : for (std::size_t i=0; i<npoin; ++i) {
1195 [ + + ]: 1699820 : if (steady) dt = m_dtp[i];
1196 [ + + ]: 10289940 : for (std::size_t c=0; c<ncomp; ++c) {
1197 : 8590120 : auto a = c*2;
1198 : 8590120 : auto b = a+1;
1199 : 8590120 : m_p(i,a) /= m_vol[i];
1200 : 8590120 : m_p(i,b) /= m_vol[i];
1201 : : // low-order solution
1202 : 8590120 : m_rhs(i,c) = m_u(i,c) - dt*m_rhs(i,c)/m_vol[i] - m_p(i,a) - m_p(i,b);
1203 : : }
1204 : : }
1205 : :
1206 : : // Allowed limits: Q+/-
1207 : :
1208 : : using std::max;
1209 : : using std::min;
1210 : :
1211 : 5330 : auto large = std::numeric_limits< tk::real >::max();
1212 [ + + ]: 1705150 : for (std::size_t i=0; i<m_q.nunk(); ++i) {
1213 [ + + ]: 10289940 : for (std::size_t c=0; c<m_q.nprop()/2; ++c) {
1214 : 8590120 : m_q(i,c*2+0) = -large;
1215 : 8590120 : m_q(i,c*2+1) = +large;
1216 : : }
1217 : : }
1218 : :
1219 : : // tetrahedron superedges
1220 [ + + ]: 984945 : for (std::size_t e=0; e<m_dsupedge[0].size()/4; ++e) {
1221 : 979615 : const auto N = m_dsupedge[0].data() + e*4;
1222 [ + + ]: 5910650 : for (std::size_t c=0; c<ncomp; ++c) {
1223 : 4931035 : auto a = c*2;
1224 : 4931035 : auto b = a+1;
1225 [ + + ]: 34517245 : for (const auto& [p,q] : tk::lpoed) {
1226 : : tk::real alwp, alwn;
1227 [ + + ]: 29586210 : if (g_cfg.get< tag::fctclip >()) {
1228 [ + - ][ + - ]: 472200 : alwp = max( m_rhs(N[p],c), m_rhs(N[q],c) );
1229 [ + - ][ + - ]: 472200 : alwn = min( m_rhs(N[p],c), m_rhs(N[q],c) );
1230 : : } else {
1231 [ + - ][ + - ]: 29114010 : alwp = max( max(m_rhs(N[p],c), m_u(N[p],c)),
1232 [ + - ][ + - ]: 29114010 : max(m_rhs(N[q],c), m_u(N[q],c)) );
1233 [ + - ][ + - ]: 29114010 : alwn = min( min(m_rhs(N[p],c), m_u(N[p],c)),
1234 [ + - ][ + - ]: 29114010 : min(m_rhs(N[q],c), m_u(N[q],c)) );
1235 : : }
1236 [ + - ][ + - ]: 29586210 : m_q(N[p],a) = max(m_q(N[p],a), alwp);
1237 [ + - ][ + - ]: 29586210 : m_q(N[p],b) = min(m_q(N[p],b), alwn);
1238 [ + - ][ + - ]: 29586210 : m_q(N[q],a) = max(m_q(N[q],a), alwp);
1239 [ + - ][ + - ]: 29586210 : m_q(N[q],b) = min(m_q(N[q],b), alwn);
1240 : : }
1241 : : }
1242 : : }
1243 : :
1244 : : // triangle superedges
1245 [ + + ]: 533300 : for (std::size_t e=0; e<m_dsupedge[1].size()/3; ++e) {
1246 : 527970 : const auto N = m_dsupedge[1].data() + e*3;
1247 [ + + ]: 3199120 : for (std::size_t c=0; c<ncomp; ++c) {
1248 : 2671150 : auto a = c*2;
1249 : 2671150 : auto b = a+1;
1250 [ + + ]: 10684600 : for (const auto& [p,q] : tk::lpoet) {
1251 : : tk::real alwp, alwn;
1252 [ + + ]: 8013450 : if (g_cfg.get< tag::fctclip >()) {
1253 [ + - ][ + - ]: 185400 : alwp = max( m_rhs(N[p],c), m_rhs(N[q],c) );
1254 [ + - ][ + - ]: 185400 : alwn = min( m_rhs(N[p],c), m_rhs(N[q],c) );
1255 : : } else {
1256 [ + - ][ + - ]: 7828050 : alwp = max( max(m_rhs(N[p],c), m_u(N[p],c)),
1257 [ + - ][ + - ]: 7828050 : max(m_rhs(N[q],c), m_u(N[q],c)) );
1258 [ + - ][ + - ]: 7828050 : alwn = min( min(m_rhs(N[p],c), m_u(N[p],c)),
1259 [ + - ][ + - ]: 7828050 : min(m_rhs(N[q],c), m_u(N[q],c)) );
1260 : : }
1261 [ + - ][ + - ]: 8013450 : m_q(N[p],a) = max(m_q(N[p],a), alwp);
1262 [ + - ][ + - ]: 8013450 : m_q(N[p],b) = min(m_q(N[p],b), alwn);
1263 [ + - ][ + - ]: 8013450 : m_q(N[q],a) = max(m_q(N[q],a), alwp);
1264 [ + - ][ + - ]: 8013450 : m_q(N[q],b) = min(m_q(N[q],b), alwn);
1265 : : }
1266 : : }
1267 : : }
1268 : :
1269 : : // edges
1270 [ + + ]: 2272950 : for (std::size_t e=0; e<m_dsupedge[2].size()/2; ++e) {
1271 : 2267620 : const auto N = m_dsupedge[2].data() + e*2;
1272 [ + + ]: 13712040 : for (std::size_t c=0; c<ncomp; ++c) {
1273 : 11444420 : auto a = c*2;
1274 : 11444420 : auto b = a+1;
1275 : : tk::real alwp, alwn;
1276 [ + + ]: 11444420 : if (g_cfg.get< tag::fctclip >()) {
1277 [ + - ][ + - ]: 214000 : alwp = max( m_rhs(N[0],c), m_rhs(N[1],c) );
1278 [ + - ][ + - ]: 214000 : alwn = min( m_rhs(N[0],c), m_rhs(N[1],c) );
1279 : : } else {
1280 [ + - ][ + - ]: 11230420 : alwp = max( max(m_rhs(N[0],c), m_u(N[0],c)),
1281 [ + - ][ + - ]: 11230420 : max(m_rhs(N[1],c), m_u(N[1],c)) );
1282 [ + - ][ + - ]: 11230420 : alwn = min( min(m_rhs(N[0],c), m_u(N[0],c)),
1283 [ + - ][ + - ]: 11230420 : min(m_rhs(N[1],c), m_u(N[1],c)) );
1284 : : }
1285 [ + - ][ + - ]: 11444420 : m_q(N[0],a) = max(m_q(N[0],a), alwp);
1286 [ + - ][ + - ]: 11444420 : m_q(N[0],b) = min(m_q(N[0],b), alwn);
1287 [ + - ][ + - ]: 11444420 : m_q(N[1],a) = max(m_q(N[1],a), alwp);
1288 [ + - ][ + - ]: 11444420 : m_q(N[1],b) = min(m_q(N[1],b), alwn);
1289 : : }
1290 : : }
1291 : :
1292 : : // Communicate allowed limits contributions
1293 [ + + ][ + + ]: 5330 : if (d->NodeCommMap().empty() or m_inactive.size() == d->NodeCommMap().size()){
[ + + ]
1294 : 145 : comalw_complete();
1295 : : } else {
1296 [ + + ]: 51830 : for (const auto& [c,n] : d->NodeCommMap()) {
1297 [ + - ][ + + ]: 46645 : if (m_inactive.count(c)) continue;
1298 : 45980 : decltype(m_qc) exp;
1299 [ + - ][ + - ]: 629110 : for (auto g : n) exp[g] = m_q[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1300 [ + - ][ + - ]: 45980 : thisProxy[c].comalw( exp );
1301 : 45980 : }
1302 : : }
1303 : 5330 : ownalw_complete();
1304 : 5330 : }
1305 : :
1306 : : void
1307 : 45980 : ZalCG::comalw( const std::unordered_map< std::size_t,
1308 : : std::vector< tk::real > >& inalw )
1309 : : // *****************************************************************************
1310 : : // Receive allowed limits contributions on chare-boundaries
1311 : : //! \param[in] inalw Partial contributions of allowed limits contributions on
1312 : : //! chare-boundary nodes. Key: global mesh node IDs, value: allowed limit
1313 : : //! contributions.
1314 : : // *****************************************************************************
1315 : : {
1316 [ + + ]: 629110 : for (const auto& [g,alw] : inalw) {
1317 [ + - ]: 583130 : auto& q = m_qc[g];
1318 [ + - ]: 583130 : q.resize( alw.size() );
1319 [ + + ]: 3543180 : for (std::size_t c=0; c<alw.size()/2; ++c) {
1320 : 2960050 : auto a = c*2;
1321 : 2960050 : auto b = a+1;
1322 : 2960050 : q[a] = std::max( q[a], alw[a] );
1323 : 2960050 : q[b] = std::min( q[b], alw[b] );
1324 : : }
1325 : : }
1326 : :
1327 : : // When we have heard from all chares we communicate with, this chare is done
1328 [ + + ]: 45980 : if (++m_nalw + m_inactive.size() == Disc()->NodeCommMap().size()) {
1329 : 5185 : m_nalw = 0;
1330 : 5185 : comalw_complete();
1331 : : }
1332 : 45980 : }
1333 : :
1334 : : void
1335 : 5330 : ZalCG::lim()
1336 : : // *****************************************************************************
1337 : : // Compute limit coefficients
1338 : : // *****************************************************************************
1339 : : {
1340 [ + - ]: 5330 : auto d = Disc();
1341 : 5330 : const auto npoin = m_u.nunk();
1342 : 5330 : const auto ncomp = m_u.nprop();
1343 : 5330 : const auto& lid = d->Lid();
1344 : :
1345 : : using std::max;
1346 : : using std::min;
1347 : :
1348 : : // Combine own and communicated contributions to allowed limits
1349 [ + + ]: 442280 : for (const auto& [g,alw] : m_qc) {
1350 [ + - ]: 436950 : auto i = tk::cref_find( lid, g );
1351 [ + + ]: 2656840 : for (std::size_t c=0; c<alw.size()/2; ++c) {
1352 : 2219890 : auto a = c*2;
1353 : 2219890 : auto b = a+1;
1354 [ + - ][ + - ]: 2219890 : m_q(i,a) = max( m_q(i,a), alw[a] );
1355 [ + - ][ + - ]: 2219890 : m_q(i,b) = min( m_q(i,b), alw[b] );
1356 : : }
1357 : : }
1358 : 5330 : tk::destroy(m_qc);
1359 : :
1360 : : // Finish computing allowed limits
1361 [ + + ]: 1705150 : for (std::size_t i=0; i<npoin; ++i) {
1362 [ + + ]: 10289940 : for (std::size_t c=0; c<ncomp; ++c) {
1363 : 8590120 : auto a = c*2;
1364 : 8590120 : auto b = a+1;
1365 [ + - ][ + - ]: 8590120 : m_q(i,a) -= m_rhs(i,c);
1366 [ + - ][ + - ]: 8590120 : m_q(i,b) -= m_rhs(i,c);
1367 : : }
1368 : : }
1369 : :
1370 : : // Limit coefficients, C
1371 : :
1372 [ + + ]: 1705150 : for (std::size_t i=0; i<npoin; ++i) {
1373 [ + + ]: 10289940 : for (std::size_t c=0; c<ncomp; ++c) {
1374 : 8590120 : auto a = c*2;
1375 : 8590120 : auto b = a+1;
1376 : 8590120 : auto eps = std::numeric_limits< tk::real >::epsilon();
1377 [ + - ][ + + ]: 8590120 : m_q(i,a) = m_p(i,a) < eps ? 0.0 : min(1.0, m_q(i,a)/m_p(i,a));
[ + - ][ + - ]
[ + - ]
1378 [ + - ][ + + ]: 8590120 : m_q(i,b) = m_p(i,b) > -eps ? 0.0 : min(1.0, m_q(i,b)/m_p(i,b));
[ + - ][ + - ]
[ + - ]
1379 : : }
1380 : : }
1381 : :
1382 : : // Limited antidiffusive contributions
1383 : :
1384 : 5330 : auto ctau = g_cfg.get< tag::fctdif >();
1385 [ + - ]: 5330 : m_a.fill( 0.0 );
1386 : :
1387 [ + - ]: 5330 : auto fctsys = g_cfg.get< tag::fctsys >();
1388 [ + + ]: 9310 : for (auto& c : fctsys) --c;
1389 : :
1390 : : #if defined(__clang__)
1391 : : #pragma clang diagnostic push
1392 : : #pragma clang diagnostic ignored "-Wvla"
1393 : : #pragma clang diagnostic ignored "-Wvla-extension"
1394 : : #elif defined(STRICT_GNUC)
1395 : : #pragma GCC diagnostic push
1396 : : #pragma GCC diagnostic ignored "-Wvla"
1397 : : #endif
1398 : :
1399 : : // tetrahedron superedges
1400 [ + + ]: 984945 : for (std::size_t e=0; e<m_dsupedge[0].size()/4; ++e) {
1401 : 979615 : const auto N = m_dsupedge[0].data() + e*4;
1402 : 979615 : const auto D = m_dsupint[0].data();
1403 : 979615 : auto C = m_dsuplim[0].data();
1404 : 979615 : std::size_t i = 0;
1405 [ + + ]: 12734995 : for (const auto& [p,q] : tk::lpoed) {
1406 : 5877690 : auto dif = D[(e*6+i)*4+3];
1407 : 5877690 : auto coef = C + (e*6+i)*ncomp;
1408 : 5877690 : tk::real aec[ncomp];
1409 [ + + ]: 35463900 : for (std::size_t c=0; c<ncomp; ++c) {
1410 [ + - ][ + - ]: 29586210 : aec[c] = -dif * ctau * (m_u(N[p],c) - m_u(N[q],c));
1411 [ + + ]: 29586210 : if (m_fctfreeze) continue;
1412 : 25835910 : auto a = c*2;
1413 : 25835910 : auto b = a+1;
1414 [ + + ][ + - ]: 25835910 : coef[c] = min( aec[c] < 0.0 ? m_q(N[p],a) : m_q(N[p],b),
[ + - ]
1415 [ + + ][ + - ]: 25835910 : aec[c] > 0.0 ? m_q(N[q],a) : m_q(N[q],b) );
[ + - ]
1416 : : }
1417 : 5877690 : tk::real cs = 1.0;
1418 [ + + ]: 18075060 : for (auto c : fctsys) cs = min( cs, coef[c] );
1419 [ + + ]: 18075060 : for (auto c : fctsys) coef[c] = cs;
1420 [ + + ]: 35463900 : for (std::size_t c=0; c<ncomp; ++c) {
1421 : 29586210 : aec[c] *= coef[c];
1422 [ + - ]: 29586210 : m_a(N[p],c) -= aec[c];
1423 [ + - ]: 29586210 : m_a(N[q],c) += aec[c];
1424 : : }
1425 : 5877690 : ++i;
1426 : 5877690 : }
1427 : : }
1428 : :
1429 : : // triangle superedges
1430 [ + + ]: 533300 : for (std::size_t e=0; e<m_dsupedge[1].size()/3; ++e) {
1431 : 527970 : const auto N = m_dsupedge[1].data() + e*3;
1432 : 527970 : const auto D = m_dsupint[1].data();
1433 : 527970 : auto C = m_dsuplim[0].data();
1434 : 527970 : std::size_t i = 0;
1435 [ + + ]: 3695790 : for (const auto& [p,q] : tk::lpoet) {
1436 : 1583910 : auto dif = D[(e*3+i)*4+3];
1437 : 1583910 : auto coef = C + (e*3+i)*ncomp;
1438 : 1583910 : tk::real aec[ncomp];
1439 [ + + ]: 9597360 : for (std::size_t c=0; c<ncomp; ++c) {
1440 [ + - ][ + - ]: 8013450 : aec[c] = -dif * ctau * (m_u(N[p],c) - m_u(N[q],c));
1441 [ + + ]: 8013450 : if (m_fctfreeze) continue;
1442 : 7255650 : auto a = c*2;
1443 : 7255650 : auto b = a+1;
1444 [ + + ][ + - ]: 7255650 : coef[c] = min( aec[c] < 0.0 ? m_q(N[p],a) : m_q(N[p],b),
[ + - ]
1445 [ + + ][ + - ]: 7255650 : aec[c] > 0.0 ? m_q(N[q],a) : m_q(N[q],b) );
[ + - ]
1446 : : }
1447 : 1583910 : tk::real cs = 1.0;
1448 [ + + ]: 5488350 : for (auto c : fctsys) cs = min( cs, coef[c] );
1449 [ + + ]: 5488350 : for (auto c : fctsys) coef[c] = cs;
1450 [ + + ]: 9597360 : for (std::size_t c=0; c<ncomp; ++c) {
1451 : 8013450 : aec[c] *= coef[c];
1452 [ + - ]: 8013450 : m_a(N[p],c) -= aec[c];
1453 [ + - ]: 8013450 : m_a(N[q],c) += aec[c];
1454 : : }
1455 : 1583910 : ++i;
1456 : 1583910 : }
1457 : : }
1458 : :
1459 : : // edges
1460 [ + + ]: 2272950 : for (std::size_t e=0; e<m_dsupedge[2].size()/2; ++e) {
1461 : 2267620 : const auto N = m_dsupedge[2].data() + e*2;
1462 : 2267620 : const auto dif = m_dsupint[2][e*4+3];
1463 : 2267620 : auto coef = m_dsuplim[2].data() + e*ncomp;
1464 : 2267620 : tk::real aec[ncomp];
1465 [ + + ]: 13712040 : for (std::size_t c=0; c<ncomp; ++c) {
1466 [ + - ][ + - ]: 11444420 : aec[c] = -dif * ctau * (m_u(N[0],c) - m_u(N[1],c));
1467 : 11444420 : auto a = c*2;
1468 : 11444420 : auto b = a+1;
1469 [ + + ]: 11444420 : if (m_fctfreeze) continue;
1470 [ + + ][ + - ]: 10176395 : coef[c] = min( aec[c] < 0.0 ? m_q(N[0],a) : m_q(N[0],b),
[ + - ]
1471 [ + + ][ + - ]: 10176395 : aec[c] > 0.0 ? m_q(N[1],a) : m_q(N[1],b) );
[ + - ]
1472 : : }
1473 : 2267620 : tk::real cs = 1.0;
1474 [ + + ]: 7348370 : for (auto c : fctsys) cs = min( cs, coef[c] );
1475 [ + + ]: 7348370 : for (auto c : fctsys) coef[c] = cs;
1476 [ + + ]: 13712040 : for (std::size_t c=0; c<ncomp; ++c) {
1477 : 11444420 : aec[c] *= coef[c];
1478 [ + - ]: 11444420 : m_a(N[0],c) -= aec[c];
1479 [ + - ]: 11444420 : m_a(N[1],c) += aec[c];
1480 : : }
1481 : 2267620 : }
1482 : :
1483 : : #if defined(__clang__)
1484 : : #pragma clang diagnostic pop
1485 : : #elif defined(STRICT_GNUC)
1486 : : #pragma GCC diagnostic pop
1487 : : #endif
1488 : :
1489 : : // Communicate limited antidiffusive contributions
1490 [ + + ][ + + ]: 5330 : if (d->NodeCommMap().empty() or m_inactive.size() == d->NodeCommMap().size()){
[ + + ]
1491 [ + - ]: 145 : comlim_complete();
1492 : : } else {
1493 [ + + ]: 51830 : for (const auto& [c,n] : d->NodeCommMap()) {
1494 [ + - ][ + + ]: 46645 : if (m_inactive.count(c)) continue;
1495 : 45980 : decltype(m_ac) exp;
1496 [ + - ][ + - ]: 629110 : for (auto g : n) exp[g] = m_a[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1497 [ + - ][ + - ]: 45980 : thisProxy[c].comlim( exp );
1498 : 45980 : }
1499 : : }
1500 [ + - ]: 5330 : ownlim_complete();
1501 : 5330 : }
1502 : :
1503 : : void
1504 : 45980 : ZalCG::comlim( const std::unordered_map< std::size_t,
1505 : : std::vector< tk::real > >& inlim )
1506 : : // *****************************************************************************
1507 : : // Receive limited antidiffusive contributions on chare-boundaries
1508 : : //! \param[in] inlim Partial contributions of limited contributions on
1509 : : //! chare-boundary nodes. Key: global mesh node IDs, value: limited
1510 : : //! contributions.
1511 : : // *****************************************************************************
1512 : : {
1513 : : using tk::operator+=;
1514 [ + - ][ + - ]: 629110 : for (const auto& [g,a] : inlim) m_ac[g] += a;
[ + + ]
1515 : :
1516 : : // When we have heard from all chares we communicate with, this chare is done
1517 [ + + ]: 45980 : if (++m_nlim + m_inactive.size() == Disc()->NodeCommMap().size()) {
1518 : 5185 : m_nlim = 0;
1519 : 5185 : comlim_complete();
1520 : : }
1521 : 45980 : }
1522 : :
1523 : : void
1524 : 6040 : ZalCG::solve()
1525 : : // *****************************************************************************
1526 : : // Advance systems of equations
1527 : : // *****************************************************************************
1528 : : {
1529 : 6040 : auto d = Disc();
1530 : 6040 : const auto npoin = m_u.nunk();
1531 : 6040 : const auto ncomp = m_u.nprop();
1532 : 6040 : const auto steady = g_cfg.get< tag::steady >();
1533 : :
1534 [ + + ]: 6040 : if (m_deactivated) {
1535 : :
1536 : 710 : m_a = m_u;
1537 : :
1538 : : } else {
1539 : :
1540 : : // Combine own and communicated contributions to limited antidiffusive
1541 : : // contributions
1542 [ + + ]: 442280 : for (const auto& [g,a] : m_ac) {
1543 [ + - ]: 436950 : auto i = tk::cref_find( d->Lid(), g );
1544 [ + - ][ + + ]: 2656840 : for (std::size_t c=0; c<a.size(); ++c) m_a(i,c) += a[c];
1545 : : }
1546 : 5330 : tk::destroy(m_ac);
1547 : :
1548 : 5330 : tk::Fields u;
1549 [ + + ]: 5330 : std::size_t cstart = m_freezeflow > 1.0 ? 5 : 0;
1550 [ + + ][ + - ]: 5330 : if (cstart) u = m_u;
1551 : :
1552 [ + - ]: 5330 : if (g_cfg.get< tag::fct >()) {
1553 : : // Apply limited antidiffusive contributions to low-order solution
1554 [ + + ]: 1705150 : for (std::size_t i=0; i<npoin; ++i) {
1555 [ + + ]: 10289940 : for (std::size_t c=0; c<ncomp; ++c) {
1556 [ + - ][ + - ]: 8590120 : m_a(i,c) = m_rhs(i,c) + m_a(i,c)/m_vol[i];
[ + - ]
1557 : : }
1558 : : }
1559 : : } else {
1560 : : // Apply rhs
1561 : 0 : auto dt = d->Dt();
1562 [ - - ]: 0 : for (std::size_t i=0; i<npoin; ++i) {
1563 [ - - ]: 0 : if (steady) dt = m_dtp[i];
1564 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c) {
1565 [ - - ][ - - ]: 0 : m_a(i,c) = m_u(i,c) - dt*m_rhs(i,c)/m_vol[i];
[ - - ]
1566 : : }
1567 : : }
1568 : : }
1569 : :
1570 : : // Configure and apply scalar source to solution (if defined)
1571 [ + - ]: 5330 : auto src = problems::PHYS_SRC();
1572 [ - + ][ - - ]: 5330 : if (src) src( d->Coord(), d->T(), m_a );
1573 : :
1574 : : // Enforce boundary conditions
1575 [ + - ]: 5330 : BC( m_a, d->T() + d->Dt() );
1576 : :
1577 : : // Explicitly zero out flow for freezeflow
1578 [ + + ]: 5330 : if (cstart) {
1579 [ + + ]: 87400 : for (std::size_t i=0; i<npoin; ++i) {
1580 [ + + ]: 518814 : for (std::size_t c=0; c<cstart; ++c) {
1581 [ + - ][ + - ]: 432345 : m_a(i,c) = u(i,c);
1582 : : }
1583 : : }
1584 : : }
1585 : :
1586 : 5330 : }
1587 : :
1588 : : // Compute diagnostics, e.g., residuals
1589 : 6040 : auto diag_iter = g_cfg.get< tag::diag_iter >();
1590 : 6040 : auto diag = m_diag.rhocompute( *d, m_a, m_u, diag_iter );
1591 : :
1592 : : // Update solution
1593 : 6040 : m_u = m_a;
1594 : 6040 : m_a.fill( 0.0 );
1595 : :
1596 : : // Increase number of iterations and physical time
1597 : 6040 : d->next();
1598 : :
1599 : : // Advance physical time for local time stepping
1600 [ + + ]: 6040 : if (steady) {
1601 : : using tk::operator+=;
1602 : 480 : m_tp += m_dtp;
1603 : : }
1604 : :
1605 : : // Evaluate residuals
1606 [ - + ][ - - ]: 6040 : if (!diag) evalres( std::vector< tk::real >( ncomp, 1.0 ) );
[ - - ]
1607 : 6040 : }
1608 : :
1609 : : void
1610 : 6040 : ZalCG::evalres( const std::vector< tk::real >& l2res )
1611 : : // *****************************************************************************
1612 : : // Evaluate residuals
1613 : : //! \param[in] l2res L2-norms of the residual for each scalar component
1614 : : //! computed across the whole problem
1615 : : // *****************************************************************************
1616 : : {
1617 : 6040 : auto d = Disc();
1618 : :
1619 [ + + ]: 6040 : if (g_cfg.get< tag::steady >()) {
1620 : 480 : const auto rc = g_cfg.get< tag::rescomp >() - 1;
1621 : 480 : d->residual( l2res[rc] );
1622 [ + + ]: 480 : if (l2res[rc] < g_cfg.get< tag::fctfreeze >()) m_fctfreeze = 1;
1623 : : }
1624 : :
1625 [ + + ]: 6040 : if (d->deactivate()) deactivate(); else refine();
1626 : 6040 : }
1627 : :
1628 : : int
1629 : 131266 : ZalCG::active( std::size_t p,
1630 : : std::size_t q,
1631 : : const std::vector< uint64_t >& sys ) const
1632 : : // *****************************************************************************
1633 : : // Decide if edge is active
1634 : : //! \param[in] p Local id of left edge-end point
1635 : : //! \param[in] q Local id of right edge-end point
1636 : : //! \param[in] sys List of components to consider
1637 : : //! \return True if active, false if not
1638 : : // *****************************************************************************
1639 : : {
1640 : 131266 : auto tol = g_cfg.get< tag::deatol >();
1641 : :
1642 : 131266 : tk::real maxdiff = 0.0;
1643 [ + + ]: 262532 : for (auto c : sys) {
1644 [ + - ][ + - ]: 131266 : auto e = std::abs( m_u(p,c) - m_u(q,c) );
1645 : 131266 : maxdiff = std::max( maxdiff, e );
1646 : : }
1647 : :
1648 : 131266 : return maxdiff > tol;
1649 : : }
1650 : :
1651 : : int
1652 : 19 : ZalCG::dea( const std::vector< uint64_t >& sys ) const
1653 : : // *****************************************************************************
1654 : : // Decide whether to deactivate this chare
1655 : : //! \param[in] sys List of components to consider
1656 : : //! \return Nonzero to deactivate, zero to keep active
1657 : : // *****************************************************************************
1658 : : {
1659 [ - + ]: 19 : if (m_toreactivate) return 0;
1660 [ - + ]: 19 : if (m_deactivated) return 1;
1661 : :
1662 : : // tetrahedron superedges
1663 [ + + ]: 7355 : for (std::size_t e=0; e<m_dsupedge[0].size()/4; ++e) {
1664 : 7337 : const auto N = m_dsupedge[0].data() + e*4;
1665 [ + + ]: 51354 : for (const auto& [p,q] : tk::lpoed) {
1666 [ + + ]: 44018 : if (active(N[p],N[q],sys)) return 0;
1667 : : }
1668 : : }
1669 : :
1670 : : // triangle superedges
1671 [ + + ]: 5175 : for (std::size_t e=0; e<m_dsupedge[1].size()/3; ++e) {
1672 : 5157 : const auto N = m_dsupedge[1].data() + e*3;
1673 [ + + ]: 20628 : for (const auto& [p,q] : tk::lpoet) {
1674 [ - + ]: 15471 : if (active(N[p],N[q],sys)) return 0;
1675 : : }
1676 : : }
1677 : :
1678 : : // edges
1679 [ + + ]: 19379 : for (std::size_t e=0; e<m_dsupedge[2].size()/2; ++e) {
1680 : 19361 : const auto N = m_dsupedge[2].data() + e*2;
1681 [ - + ]: 19361 : if (active(N[0],N[1],sys)) return 0;
1682 : : }
1683 : :
1684 : 18 : return 1;
1685 : : }
1686 : :
1687 : : void
1688 : 48 : ZalCG::rea( const std::vector< uint64_t >& sys, std::unordered_set< int >& req )
1689 : : const
1690 : : // *****************************************************************************
1691 : : // Decide whether to reactivate any neighbor chare
1692 : : //! \param[in] sys List of components to consider
1693 : : //! \param[in,out] req Set of neighbor chares to reactivate
1694 : : // *****************************************************************************
1695 : : {
1696 [ + + ]: 401 : for (const auto& [c,edges] : m_chbndedge) {
1697 [ + - ][ + + ]: 353 : if (not m_inactive.count(c)) continue;
1698 [ + + ]: 52547 : for (const auto& [e,sint] : edges) {
1699 [ + - ][ + + ]: 52416 : if (active(e[0],e[1],sys)) {
1700 [ + - ]: 6 : req.insert(c);
1701 : 6 : break;
1702 : : }
1703 : : }
1704 : : }
1705 : 48 : }
1706 : :
1707 : : void
1708 : 48 : ZalCG::huldif()
1709 : : // *****************************************************************************
1710 : : // Apply diffusion on active hull
1711 : : // *****************************************************************************
1712 : : {
1713 : 48 : const auto eps = std::numeric_limits< tk::real >::epsilon();
1714 : 48 : const auto dif = g_cfg.get< tag::deadif >();
1715 [ - + ]: 48 : if (std::abs(dif) < eps) return;
1716 : :
1717 : 48 : const auto npoin = m_u.nunk();
1718 : 48 : const auto ncomp = m_u.nprop();
1719 : :
1720 : 48 : m_rhs.fill( 0.0 );
1721 [ + + ]: 401 : for (const auto& [ch,edges] : m_chbndedge) {
1722 [ + - ][ + + ]: 353 : if (not m_inactive.count(ch)) continue;
1723 [ + + ]: 56673 : for (const auto& [e,sint] : edges) {
1724 : 56536 : auto p = e[0];
1725 : 56536 : auto q = e[1];
1726 : 56536 : auto fw = dif * sint;
1727 [ + + ]: 339216 : for (std::size_t c=0; c<ncomp; ++c) {
1728 [ + - ][ + - ]: 282680 : auto f = fw*(m_u(p,c) - m_u(q,c));
1729 [ + - ]: 282680 : m_rhs(p,c) -= f;
1730 [ + - ]: 282680 : m_rhs(q,c) += f;
1731 : : }
1732 : : }
1733 : : }
1734 : :
1735 [ + + ]: 37506 : for (std::size_t i=0; i<npoin; ++i) {
1736 [ + + ]: 224748 : for (std::size_t c=0; c<ncomp; ++c) {
1737 : 187290 : m_u(i,c) += m_rhs(i,c)/m_vol[i];
1738 : : }
1739 : : }
1740 : : }
1741 : :
1742 : : void
1743 : 190 : ZalCG::deactivate()
1744 : : // *****************************************************************************
1745 : : // Deactivate regions
1746 : : // *****************************************************************************
1747 : : {
1748 [ + - ]: 190 : auto d = Disc();
1749 : :
1750 : 190 : std::unordered_set< int > reactivate;
1751 : :
1752 [ + + ]: 190 : if (not m_deactivated) {
1753 : : // Apply diffusion on active hull
1754 [ + - ]: 48 : huldif();
1755 : : // Query scalar components to deactivate based on
1756 [ + - ]: 48 : auto sys = g_cfg.get< tag::deasys >();
1757 [ + + ]: 96 : for (auto& c : sys) --c;
1758 : : // Decide whether to deactivate this chare
1759 [ + - ][ + + ]: 48 : if (d->deastart() and dea(sys)) m_todeactivate = 1;
[ + - ][ + + ]
[ + + ]
1760 : : // Decide whether to reactivate any neighbor chare
1761 [ + - ]: 48 : rea( sys, reactivate );
1762 : 48 : }
1763 : :
1764 : : // Communicate deactivation requests
1765 [ + + ]: 1550 : for (const auto& [c,n] : d->NodeCommMap()) {
1766 [ + - ][ + - ]: 1360 : thisProxy[c].comdea( reactivate.count(c) );
[ + - ]
1767 : : }
1768 [ + - ]: 190 : owndea_complete();
1769 : 190 : }
1770 : :
1771 : : void
1772 : 1360 : ZalCG::comdea( std::size_t reactivate )
1773 : : // *****************************************************************************
1774 : : // Communicate deactivation desires
1775 : : //! \param[in] reactivate 1 if sender wants to reactivate this chare, 0 if not
1776 : : // *****************************************************************************
1777 : : {
1778 : 1360 : auto d = Disc();
1779 : :
1780 [ + + ][ + + ]: 1360 : if (m_deactivated and reactivate) m_toreactivate = 1;
1781 : :
1782 : : // Communicate activation status to neighbors
1783 [ + + ]: 1360 : if (++m_ndea == d->NodeCommMap().size()) {
1784 : 190 : m_ndea = 0;
1785 : 190 : comdea_complete();
1786 : : }
1787 : 1360 : }
1788 : :
1789 : : void
1790 : 190 : ZalCG::activate()
1791 : : // *****************************************************************************
1792 : : // Compute deactivation status
1793 : : // *****************************************************************************
1794 : : {
1795 : 190 : auto d = Disc();
1796 : :
1797 [ + + ]: 190 : if (m_todeactivate) m_deactivated = 1;
1798 [ + + ]: 190 : if (m_toreactivate) m_deactivated = 0;
1799 : :
1800 : : // Communicate deactivation status
1801 [ + + ]: 1550 : for (const auto& [c,n] : d->NodeCommMap()) {
1802 [ + - ][ + - ]: 1360 : thisProxy[c].comact( thisIndex, m_deactivated );
1803 : : }
1804 : 190 : ownact_complete();
1805 : 190 : }
1806 : :
1807 : : void
1808 : 1360 : ZalCG::comact( int ch, int deactivated )
1809 : : // *****************************************************************************
1810 : : // Communicate deactivation status
1811 : : //! \param[in] ch Sender chare id
1812 : : //! \param[in] deactivated 1 if sender was deactivated, 0 if not
1813 : : // *****************************************************************************
1814 : : {
1815 : 1360 : auto d = Disc();
1816 : :
1817 [ + + ]: 1360 : if (deactivated) m_inactive.insert(ch); else m_inactive.erase(ch);
1818 : :
1819 : : // When we have heard from all chares we communicate with, this chare is done
1820 [ + + ]: 1360 : if (++m_nact == d->NodeCommMap().size()) {
1821 : : // Aggregate deactivation status to every chare
1822 [ + - ]: 190 : contribute( sizeof(int), &m_deactivated, CkReduction::sum_int,
1823 [ + - ][ + - ]: 380 : CkCallback(CkReductionTarget(ZalCG,deastat), thisProxy) );
1824 : 190 : m_todeactivate = 0;
1825 : 190 : m_toreactivate = 0;
1826 : 190 : m_nact = 0;
1827 : 190 : comact_complete();
1828 : : }
1829 : 1360 : }
1830 : :
1831 : : void
1832 : 190 : ZalCG::deastat( int dea )
1833 : : // *****************************************************************************
1834 : : // Receive aggregate deactivation status
1835 : : //! \param[in] dea Sum of deactived chares
1836 : : // *****************************************************************************
1837 : : {
1838 : 190 : auto d = Disc();
1839 : :
1840 : : // Disallow deactivating all chares
1841 [ - + ]: 190 : if (dea == d->nchare()) {
1842 : 0 : dea = 0;
1843 : 0 : m_deactivated = 0;
1844 : 0 : m_inactive.clear();
1845 : : }
1846 : :
1847 : : // Report number of deactivated chares for diagnostics
1848 [ + + ]: 190 : if (thisIndex == 0) d->deactivated( dea );
1849 : :
1850 : 190 : deastat_complete();
1851 : 190 : }
1852 : :
1853 : : void
1854 : 6040 : ZalCG::refine()
1855 : : // *****************************************************************************
1856 : : // Optionally refine/derefine mesh
1857 : : // *****************************************************************************
1858 : : {
1859 : 6040 : auto d = Disc();
1860 : :
1861 : : // See if this is the last time step
1862 [ + + ]: 6040 : if (d->finished()) m_finished = 1;
1863 : :
1864 : 6040 : const auto& ht = g_cfg.get< tag::href >();
1865 : 6040 : auto dtref = ht.get< tag::dt >();
1866 : 6040 : auto dtfreq = ht.get< tag::dtfreq >();
1867 : :
1868 : : // if t>0 refinement enabled and we hit the frequency
1869 [ - + ][ - - ]: 6040 : if (dtref && !(d->It() % dtfreq)) { // refine
[ - + ]
1870 : :
1871 : 0 : d->refined() = 1;
1872 : 0 : d->startvol();
1873 : 0 : d->Ref()->dtref( m_bface, m_bnode, m_triinpoel );
1874 : :
1875 : : // Activate SDAG waits for re-computing the integrals
1876 [ - - ][ - - ]: 0 : thisProxy[ thisIndex ].wait4int();
1877 : :
1878 : : } else { // do not refine
1879 : :
1880 : 6040 : d->refined() = 0;
1881 : 6040 : feop_complete();
1882 : 6040 : resize_complete();
1883 : :
1884 : : }
1885 : 6040 : }
1886 : :
1887 : : void
1888 : 0 : ZalCG::resizePostAMR(
1889 : : const std::vector< std::size_t >& /*ginpoel*/,
1890 : : const tk::UnsMesh::Chunk& chunk,
1891 : : const tk::UnsMesh::Coords& coord,
1892 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
1893 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
1894 : : const std::set< std::size_t >& removedNodes,
1895 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
1896 : : nodeCommMap,
1897 : : const std::map< int, std::vector< std::size_t > >& bface,
1898 : : const std::map< int, std::vector< std::size_t > >& bnode,
1899 : : const std::vector< std::size_t >& triinpoel )
1900 : : // *****************************************************************************
1901 : : // Receive new mesh from Refiner
1902 : : //! \param[in] ginpoel Mesh connectivity with global node ids
1903 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
1904 : : //! \param[in] coord New mesh node coordinates
1905 : : //! \param[in] addedNodes Newly added mesh nodes and their parents (local ids)
1906 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1907 : : //! \param[in] removedNodes Newly removed mesh node local ids
1908 : : //! \param[in] nodeCommMap New node communication map
1909 : : //! \param[in] bface Boundary-faces mapped to side set ids
1910 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
1911 : : //! \param[in] triinpoel Boundary-face connectivity
1912 : : // *****************************************************************************
1913 : : {
1914 [ - - ]: 0 : auto d = Disc();
1915 : 0 : auto ncomp = m_u.nprop();
1916 : :
1917 : 0 : d->Itf() = 0; // Zero field output iteration count if AMR
1918 : 0 : ++d->Itr(); // Increase number of iterations with a change in the mesh
1919 : :
1920 : : // Resize mesh data structures after mesh refinement
1921 [ - - ]: 0 : d->resizePostAMR( chunk, coord, nodeCommMap, removedNodes );
1922 : :
1923 [ - - ][ - - ]: 0 : Assert(coord[0].size() == m_u.nunk()-removedNodes.size()+addedNodes.size(),
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
1924 : : "Incorrect vector length post-AMR: expected length after resizing = " +
1925 : : std::to_string(coord[0].size()) + ", actual unknown vector length = " +
1926 : : std::to_string(m_u.nunk()-removedNodes.size()+addedNodes.size()));
1927 : :
1928 : : // Remove newly removed nodes from solution vectors
1929 [ - - ]: 0 : m_u.rm( removedNodes );
1930 [ - - ]: 0 : m_rhs.rm( removedNodes );
1931 : :
1932 : : // Resize auxiliary solution vectors
1933 : 0 : auto npoin = coord[0].size();
1934 [ - - ]: 0 : m_u.resize( npoin );
1935 [ - - ]: 0 : m_rhs.resize( npoin );
1936 : :
1937 : : // Update solution on new mesh
1938 [ - - ]: 0 : for (const auto& n : addedNodes)
1939 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c) {
1940 [ - - ][ - - ]: 0 : Assert(n.first < m_u.nunk(), "Added node index out of bounds post-AMR");
[ - - ][ - - ]
1941 [ - - ][ - - ]: 0 : Assert(n.second[0] < m_u.nunk() && n.second[1] < m_u.nunk(),
[ - - ][ - - ]
[ - - ]
1942 : : "Indices of parent-edge nodes out of bounds post-AMR");
1943 [ - - ][ - - ]: 0 : m_u(n.first,c) = (m_u(n.second[0],c) + m_u(n.second[1],c))/2.0;
[ - - ]
1944 : : }
1945 : :
1946 : : // Update physical-boundary node-, face-, and element lists
1947 [ - - ]: 0 : m_bnode = bnode;
1948 [ - - ]: 0 : m_bface = bface;
1949 [ - - ]: 0 : m_triinpoel = tk::remap( triinpoel, d->Lid() );
1950 : :
1951 : 0 : auto meshid = d->MeshId();
1952 [ - - ]: 0 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1953 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,resized), d->Tr()) );
1954 : 0 : }
1955 : :
1956 : : void
1957 : 1181 : ZalCG::writeFields( CkCallback cb )
1958 : : // *****************************************************************************
1959 : : // Output mesh-based fields to file
1960 : : //! \param[in] cb Function to continue with after the write
1961 : : // *****************************************************************************
1962 : : {
1963 [ + + ][ + - ]: 1181 : if (g_cfg.get< tag::benchmark >()) { cb.send(); return; }
1964 : :
1965 [ + - ]: 597 : auto d = Disc();
1966 : 597 : auto ncomp = m_u.nprop();
1967 : 597 : auto steady = g_cfg.get< tag::steady >();
1968 : :
1969 : : // Field output
1970 : :
1971 : : std::vector< std::string > nodefieldnames
1972 [ + - ][ + + ]: 4179 : {"density", "velocityx", "velocityy", "velocityz", "energy", "pressure"};
[ - - ]
1973 [ + - ][ + - ]: 72 : if (steady) nodefieldnames.push_back( "mach" );
1974 : :
1975 : : using tk::operator/=;
1976 [ + - ]: 597 : auto r = m_u.extract(0);
1977 [ + - ][ + - ]: 597 : auto u = m_u.extract(1); u /= r;
1978 [ + - ][ + - ]: 597 : auto v = m_u.extract(2); v /= r;
1979 [ + - ][ + - ]: 597 : auto w = m_u.extract(3); w /= r;
1980 [ + - ][ + - ]: 597 : auto e = m_u.extract(4); e /= r;
1981 [ + - ]: 597 : std::vector< tk::real > pr( m_u.nunk() ), ma;
1982 [ + + ][ + - ]: 597 : if (steady) ma.resize( m_u.nunk() );
1983 [ + + ]: 437228 : for (std::size_t i=0; i<pr.size(); ++i) {
1984 : 436631 : auto vv = u[i]*u[i] + v[i]*v[i] + w[i]*w[i];
1985 : 436631 : pr[i] = eos::pressure( r[i]*(e[i] - 0.5*vv) );
1986 [ + + ]: 436631 : if (steady) ma[i] = std::sqrt(vv) / eos::soundspeed( r[i], pr[i] );
1987 : : }
1988 : :
1989 : : std::vector< std::vector< tk::real > > nodefields{
1990 : 2985 : std::move(r), std::move(u), std::move(v), std::move(w), std::move(e),
1991 [ + - ][ + + ]: 4179 : std::move(pr) };
[ - - ]
1992 [ + + ][ + - ]: 597 : if (steady) nodefields.push_back( std::move(ma) );
1993 : :
1994 [ + + ]: 695 : for (std::size_t c=0; c<ncomp-5; ++c) {
1995 [ + - ][ + - ]: 98 : nodefieldnames.push_back( "c" + std::to_string(c) );
[ + - ]
1996 [ + - ][ + - ]: 98 : nodefields.push_back( m_u.extract(5+c) );
1997 : : }
1998 : :
1999 : : // query function to evaluate analytic solution (if defined)
2000 [ + - ]: 597 : auto sol = problems::SOL();
2001 : :
2002 [ + + ]: 597 : if (sol) {
2003 : 153 : const auto& coord = d->Coord();
2004 : 153 : const auto& x = coord[0];
2005 : 153 : const auto& y = coord[1];
2006 : 153 : const auto& z = coord[2];
2007 [ + - ]: 153 : auto an = m_u;
2008 [ + - ]: 153 : std::vector< tk::real > ap( m_u.nunk() );
2009 [ + + ]: 15382 : for (std::size_t i=0; i<an.nunk(); ++i) {
2010 [ + - ]: 15229 : auto s = sol( x[i], y[i], z[i], d->T(), /*meshid=*/0 );
2011 : 15229 : s[1] /= s[0];
2012 : 15229 : s[2] /= s[0];
2013 : 15229 : s[3] /= s[0];
2014 : 15229 : s[4] /= s[0];
2015 [ + - ][ + + ]: 100476 : for (std::size_t c=0; c<s.size(); ++c) an(i,c) = s[c];
2016 : 15229 : s[4] -= 0.5*(s[1]*s[1] + s[2]*s[2] + s[3]*s[3]);
2017 : 15229 : ap[i] = eos::pressure( s[0]*s[4] );
2018 : 15229 : }
2019 [ + + ]: 918 : for (std::size_t c=0; c<5; ++c) {
2020 [ + - ][ + - ]: 765 : nodefieldnames.push_back( nodefieldnames[c] + "_analytic" );
2021 [ + - ][ + - ]: 765 : nodefields.push_back( an.extract(c) );
2022 : : }
2023 [ + - ][ + - ]: 153 : nodefieldnames.push_back( nodefieldnames[5] + "_analytic" );
2024 [ + - ]: 153 : nodefields.push_back( std::move(ap) );
2025 [ + + ]: 251 : for (std::size_t c=0; c<ncomp-5; ++c) {
2026 [ + - ][ + - ]: 98 : nodefieldnames.push_back( nodefieldnames[6+c] + "_analytic" );
2027 [ + - ][ + - ]: 98 : nodefields.push_back( an.extract(5+c) );
2028 : : }
2029 : 153 : }
2030 : :
2031 : 597 : std::vector< tk::real > bnded, hull, dea;
2032 [ + + ]: 597 : if (g_cfg.get< tag::deactivate >()) {
2033 [ + - ][ + - ]: 209 : nodefieldnames.push_back( "chbnded" );
2034 [ + - ]: 209 : bnded.resize( m_u.nunk(), 0.0 );
2035 [ + + ]: 1705 : for (const auto& [ch,edges] : m_chbndedge) {
2036 [ + + ]: 794233 : for (const auto& [ed,sint] : edges) {
2037 : 792737 : bnded[ed[0]] = bnded[ed[1]] = 1.0;
2038 : : }
2039 : : }
2040 [ + - ]: 209 : nodefields.push_back( bnded );
2041 : :
2042 [ + - ][ + - ]: 209 : nodefieldnames.push_back( "activehull" );
2043 [ + - ]: 209 : hull.resize( m_u.nunk(), 0.0 );
2044 [ + + ]: 1705 : for (const auto& [ch,edges] : m_chbndedge) {
2045 [ + - ][ + + ]: 1496 : if (m_inactive.count(ch)) {
2046 [ + + ]: 590585 : for (const auto& [ed,sint] : edges) {
2047 : 589491 : hull[ed[0]] = hull[ed[1]] = 1.0;
2048 : : }
2049 : : }
2050 : : }
2051 [ + - ]: 209 : nodefields.push_back( hull );
2052 : :
2053 [ + - ][ + - ]: 209 : nodefieldnames.push_back( "deactivated" );
2054 [ + - ]: 209 : dea.resize( m_u.nunk(), static_cast<tk::real>(m_deactivated) );
2055 [ + - ]: 209 : nodefields.push_back( dea );
2056 : : }
2057 : :
2058 [ - + ][ - - ]: 597 : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
[ - - ][ - - ]
2059 : :
2060 : : // Surface output
2061 : :
2062 : 597 : std::vector< std::string > nodesurfnames;
2063 : 597 : std::vector< std::vector< tk::real > > nodesurfs;
2064 : :
2065 : 597 : const auto& f = g_cfg.get< tag::fieldout, tag::sidesets >();
2066 : :
2067 [ + + ]: 597 : if (!f.empty()) {
2068 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "density" );
2069 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "velocityx" );
2070 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "velocityy" );
2071 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "velocityz" );
2072 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "energy" );
2073 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "pressure" );
2074 : :
2075 [ - + ]: 19 : for (std::size_t c=0; c<ncomp-5; ++c) {
2076 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "c" + std::to_string(c) );
[ - - ]
2077 : : }
2078 : :
2079 [ - + ]: 19 : if (g_cfg.get< tag::deactivate >()) {
2080 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "chbnded" );
2081 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "activehull" );
2082 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "deactivated" );
2083 : : }
2084 : :
2085 [ - + ][ - - ]: 19 : if (steady) nodesurfnames.push_back( "mach" );
[ - - ]
2086 : :
2087 [ + - ]: 19 : auto bnode = tk::bfacenodes( m_bface, m_triinpoel );
2088 [ + - ]: 19 : std::set< int > outsets( begin(f), end(f) );
2089 [ + + ]: 71 : for (auto sideset : outsets) {
2090 [ + - ]: 52 : auto b = bnode.find(sideset);
2091 [ + + ]: 52 : if (b == end(bnode)) continue;
2092 : 46 : const auto& nodes = b->second;
2093 : 46 : auto i = nodesurfs.size();
2094 : 46 : auto ns = ncomp + 1;
2095 [ - + ]: 46 : if (g_cfg.get< tag::deactivate >()) ns += 3;
2096 [ - + ]: 46 : if (steady) ++ns;
2097 [ + - ]: 46 : nodesurfs.insert( end(nodesurfs), ns,
2098 [ + - ]: 92 : std::vector< tk::real >( nodes.size() ) );
2099 : 46 : std::size_t j = 0;
2100 [ + + ]: 4876 : for (auto n : nodes) {
2101 [ + - ]: 4830 : const auto s = m_u[n];
2102 : 4830 : std::size_t p = 0;
2103 : 4830 : nodesurfs[i+(p++)][j] = s[0];
2104 : 4830 : nodesurfs[i+(p++)][j] = s[1]/s[0];
2105 : 4830 : nodesurfs[i+(p++)][j] = s[2]/s[0];
2106 : 4830 : nodesurfs[i+(p++)][j] = s[3]/s[0];
2107 : 4830 : nodesurfs[i+(p++)][j] = s[4]/s[0];
2108 : 4830 : auto vv = (s[1]*s[1] + s[2]*s[2] + s[3]*s[3])/s[0]/s[0];
2109 : 4830 : auto ei = s[4]/s[0] - 0.5*vv;
2110 : 4830 : auto sp = eos::pressure( s[0]*ei );
2111 : 4830 : nodesurfs[i+(p++)][j] = sp;
2112 [ - + ]: 4830 : for (std::size_t c=0; c<ncomp-5; ++c) nodesurfs[i+(p++)+c][j] = s[5+c];
2113 [ - + ]: 4830 : if (g_cfg.get< tag::deactivate >()) {
2114 : 0 : nodesurfs[i+(p++)][j] = bnded[n];
2115 : 0 : nodesurfs[i+(p++)][j] = hull[n];
2116 : 0 : nodesurfs[i+(p++)][j] = dea[n];
2117 : : }
2118 [ - + ]: 4830 : if (steady) {
2119 : 0 : nodesurfs[i+(p++)][j] = std::sqrt(vv) / eos::soundspeed( s[0], sp );
2120 : : }
2121 : 4830 : ++j;
2122 : 4830 : }
2123 : : }
2124 : 19 : }
2125 : :
2126 : : // Send mesh and fields data (solution dump) for output to file
2127 [ + - ][ + - ]: 1194 : d->write( d->Inpoel(), d->Coord(), m_bface, tk::remap(m_bnode,d->Lid()),
2128 : 597 : m_triinpoel, {}, nodefieldnames, {}, nodesurfnames,
2129 : : {}, nodefields, {}, nodesurfs, cb );
2130 [ + - ][ + - ]: 2388 : }
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ][ - - ]
[ - - ][ - - ]
[ - - ]
2131 : :
2132 : : void
2133 : 6040 : ZalCG::out()
2134 : : // *****************************************************************************
2135 : : // Output mesh field data
2136 : : // *****************************************************************************
2137 : : {
2138 : 6040 : auto d = Disc();
2139 : :
2140 : : // Time history
2141 [ + + ][ + + ]: 6040 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ][ + + ]
2142 : 214 : auto ncomp = m_u.nprop();
2143 : 214 : const auto& inpoel = d->Inpoel();
2144 [ + - ]: 214 : std::vector< std::vector< tk::real > > hist( d->Hist().size() );
2145 : 214 : std::size_t j = 0;
2146 [ + + ]: 282 : for (const auto& p : d->Hist()) {
2147 : 68 : auto e = p.get< tag::elem >(); // host element id
2148 : 68 : const auto& n = p.get< tag::fn >(); // shapefunctions evaluated at point
2149 [ + - ]: 68 : hist[j].resize( ncomp+1, 0.0 );
2150 [ + + ]: 340 : for (std::size_t i=0; i<4; ++i) {
2151 [ + - ]: 272 : const auto u = m_u[ inpoel[e*4+i] ];
2152 : 272 : hist[j][0] += n[i] * u[0];
2153 : 272 : hist[j][1] += n[i] * u[1]/u[0];
2154 : 272 : hist[j][2] += n[i] * u[2]/u[0];
2155 : 272 : hist[j][3] += n[i] * u[3]/u[0];
2156 : 272 : hist[j][4] += n[i] * u[4]/u[0];
2157 : 272 : auto ei = u[4]/u[0] - 0.5*(u[1]*u[1] + u[2]*u[2] + u[3]*u[3])/u[0]/u[0];
2158 : 272 : hist[j][5] += n[i] * eos::pressure( u[0]*ei );
2159 [ - + ]: 272 : for (std::size_t c=5; c<ncomp; ++c) hist[j][c+1] += n[i] * u[c];
2160 : 272 : }
2161 : 68 : ++j;
2162 : : }
2163 [ + - ]: 214 : d->history( std::move(hist) );
2164 : 214 : }
2165 : :
2166 : : // Field data
2167 [ + + ][ + + ]: 6040 : if (d->fielditer() or d->fieldtime() or d->fieldrange() or m_finished) {
[ + - ][ + + ]
[ + + ]
2168 [ + - ][ + - ]: 746 : writeFields( CkCallback(CkIndex_ZalCG::integrals(), thisProxy[thisIndex]) );
[ + - ][ + - ]
2169 : : } else {
2170 : 5294 : integrals();
2171 : : }
2172 : 6040 : }
2173 : :
2174 : : void
2175 : 6040 : ZalCG::integrals()
2176 : : // *****************************************************************************
2177 : : // Compute integral quantities for output
2178 : : // *****************************************************************************
2179 : : {
2180 : 6040 : auto d = Disc();
2181 : :
2182 [ + - ][ + - ]: 6040 : if (d->integiter() or d->integtime() or d->integrange()) {
[ - + ][ - + ]
2183 : :
2184 : : using namespace integrals;
2185 [ - - ]: 0 : std::vector< std::map< int, tk::real > > ints( NUMINT );
2186 : :
2187 : : // Prepend integral vector with metadata on the current time step:
2188 : : // current iteration count, current physical time, time step size
2189 [ - - ]: 0 : ints[ ITER ][ 0 ] = static_cast< tk::real >( d->It() );
2190 [ - - ]: 0 : ints[ TIME ][ 0 ] = d->T();
2191 [ - - ]: 0 : ints[ DT ][ 0 ] = d->Dt();
2192 : :
2193 : : // Compute integrals requested for surfaces requested
2194 : 0 : const auto& reqv = g_cfg.get< tag::integout, tag::integrals >();
2195 [ - - ]: 0 : std::unordered_set< std::string > req( begin(reqv), end(reqv) );
2196 [ - - ][ - - ]: 0 : if (req.count("mass_flow_rate")) {
[ - - ]
2197 [ - - ]: 0 : for (const auto& [s,sint] : m_surfint) {
2198 [ - - ]: 0 : auto& mfr = ints[ MASS_FLOW_RATE ][ s ];
2199 : 0 : const auto& nodes = sint.first;
2200 : 0 : const auto& ndA = sint.second;
2201 : 0 : auto n = ndA.data();
2202 [ - - ]: 0 : for (auto p : nodes) {
2203 [ - - ][ - - ]: 0 : mfr += n[0]*m_u(p,1) + n[1]*m_u(p,2) + n[2]*m_u(p,3);
[ - - ]
2204 : 0 : n += 3;
2205 : : }
2206 : : }
2207 : : }
2208 : :
2209 [ - - ]: 0 : auto stream = serialize( d->MeshId(), ints );
2210 [ - - ]: 0 : d->contribute( stream.first, stream.second.get(), IntegralsMerger,
2211 [ - - ][ - - ]: 0 : CkCallback(CkIndex_Transporter::integrals(nullptr), d->Tr()) );
2212 : :
2213 : 0 : } else {
2214 : :
2215 : 6040 : step();
2216 : :
2217 : : }
2218 : 6040 : }
2219 : :
2220 : : void
2221 : 5605 : ZalCG::evalLB( int nrestart )
2222 : : // *****************************************************************************
2223 : : // Evaluate whether to do load balancing
2224 : : //! \param[in] nrestart Number of times restarted
2225 : : // *****************************************************************************
2226 : : {
2227 : 5605 : auto d = Disc();
2228 : :
2229 : : // Detect if just returned from a checkpoint and if so, zero timers and
2230 : : // finished flag
2231 [ + + ]: 5605 : if (d->restarted( nrestart )) m_finished = 0;
2232 : :
2233 : : // Load balancing if user frequency is reached or after the second time-step
2234 [ + + ]: 5605 : if (d->lb()) {
2235 : :
2236 : 3614 : AtSync();
2237 [ - + ]: 3614 : if (g_cfg.get< tag::nonblocking >()) dt();
2238 : :
2239 : : } else {
2240 : :
2241 : 1991 : dt();
2242 : :
2243 : : }
2244 : 5605 : }
2245 : :
2246 : : void
2247 : 5600 : ZalCG::evalRestart()
2248 : : // *****************************************************************************
2249 : : // Evaluate whether to save checkpoint/restart
2250 : : // *****************************************************************************
2251 : : {
2252 : 5600 : auto d = Disc();
2253 : :
2254 : 5600 : const auto rsfreq = g_cfg.get< tag::rsfreq >();
2255 : 5600 : const auto benchmark = g_cfg.get< tag::benchmark >();
2256 : :
2257 [ + + ][ - + ]: 5600 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
[ - + ]
2258 : :
2259 [ - - ]: 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
2260 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
2261 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
2262 : :
2263 : 0 : } else {
2264 : :
2265 : 5600 : evalLB( /* nrestart = */ -1 );
2266 : :
2267 : : }
2268 : 5600 : }
2269 : :
2270 : : void
2271 : 6040 : ZalCG::step()
2272 : : // *****************************************************************************
2273 : : // Evaluate whether to continue with next time step
2274 : : // *****************************************************************************
2275 : : {
2276 : 6040 : auto d = Disc();
2277 : :
2278 : : // Output one-liner status report to screen
2279 [ + + ]: 6040 : if (thisIndex == 0) d->status();
2280 : :
2281 [ + + ]: 6040 : if (not m_finished) {
2282 : :
2283 : 5600 : evalRestart();
2284 : :
2285 : : } else {
2286 : :
2287 : 440 : auto meshid = d->MeshId();
2288 [ + - ]: 440 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
2289 [ + - ][ + - ]: 880 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
2290 : :
2291 : : }
2292 : 6040 : }
2293 : :
2294 : : #include "NoWarning/zalcg.def.h"
|
