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