Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/ZalCG.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.,
7 : : 2022-2025 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief ZalCG: Taylor-Galerkin, FCT, edge-based continuous Galerkin
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "XystBuildConfig.hpp"
14 : : #include "ZalCG.hpp"
15 : : #include "Vector.hpp"
16 : : #include "Reader.hpp"
17 : : #include "ContainerUtil.hpp"
18 : : #include "UnsMesh.hpp"
19 : : #include "ExodusIIMeshWriter.hpp"
20 : : #include "InciterConfig.hpp"
21 : : #include "DerivedData.hpp"
22 : : #include "Discretization.hpp"
23 : : #include "DiagReducer.hpp"
24 : : #include "IntegralReducer.hpp"
25 : : #include "Integrals.hpp"
26 : : #include "Refiner.hpp"
27 : : #include "Reorder.hpp"
28 : : #include "Around.hpp"
29 : : #include "Zalesak.hpp"
30 : : #include "Problems.hpp"
31 : : #include "EOS.hpp"
32 : : #include "BC.hpp"
33 : :
34 : : namespace inciter {
35 : :
36 : : extern ctr::Config g_cfg;
37 : :
38 : : static CkReduction::reducerType IntegralsMerger;
39 : :
40 : : } // inciter::
41 : :
42 : : using inciter::g_cfg;
43 : : using inciter::ZalCG;
44 : :
45 : 435 : ZalCG::ZalCG( const CProxy_Discretization& disc,
46 : : const std::map< int, std::vector< std::size_t > >& bface,
47 : : const std::map< int, std::vector< std::size_t > >& bnode,
48 : 435 : const std::vector< std::size_t >& triinpoel ) :
49 : : 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 [ + - ]: 295 : auto& n = sym[ k->first ];
210 [ + - ][ + + ]: 35706 : 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 [ + + ]: 730 : 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 [ + + ]: 4578 : for (const auto& [c,nodes] : d->NodeCommMap()) {
268 : : decltype(m_bnorm) exp;
269 [ + + ]: 45132 : for (auto i : nodes) {
270 [ + + ]: 113900 : for (const auto& [s,b] : m_bnorm) {
271 : : auto k = b.find(i);
272 [ + + ]: 86419 : if (k != end(b)) exp[s][i] = k->second;
273 : : }
274 : : }
275 [ + - ][ + - ]: 8308 : 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 : 4154 : 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 [ + + ]: 7479 : for (const auto& [s,b] : inbnd) {
410 : : auto& bndnorm = m_bnormc[s];
411 [ + + ]: 16822 : for (const auto& [p,n] : b) {
412 : : auto& norm = bndnorm[p];
413 : 13497 : norm[0] += n[0];
414 : 13497 : norm[1] += n[1];
415 : 13497 : norm[2] += n[2];
416 : 13497 : norm[3] += n[3];
417 : : }
418 : : }
419 : :
420 [ + + ]: 4154 : if (++m_nnorm == Disc()->NodeCommMap().size()) {
421 : 424 : m_nnorm = 0;
422 : 424 : comnorm_complete();
423 : : }
424 : 4154 : }
425 : :
426 : : void
427 : 784 : 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 : 784 : NodeDiagnostics::registerReducers();
438 : 784 : IntegralsMerger = CkReduction::addReducer( integrals::mergeIntegrals );
439 : 784 : }
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 [ + + ]: 11608 : for (const auto& [g,n] : b) {
512 : : auto& norm = bndnorm[g];
513 : 10677 : norm[0] += n[0];
514 : 10677 : norm[1] += n[1];
515 : 10677 : norm[2] += n[2];
516 : 10677 : 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 [ + + ]: 39179 : for (auto& [g,n] : b) {
524 : 38180 : n[0] /= n[3];
525 : 38180 : n[1] /= n[3];
526 : 38180 : 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 [ + + ]: 39179 : for (auto&& [g,n] : b) {
538 : 38180 : 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 [ + + ]: 34264 : for (const auto& [g,b] : m_bndpoinint) {
557 : 33829 : m_bpoin[i] = tk::cref_find( lid, g );
558 : 33829 : m_bpint[i*3+0] = b[0];
559 : 33829 : m_bpint[i*3+1] = b[1];
560 : 33829 : m_bpint[i*3+2] = b[2];
561 : 33829 : ++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 [ + + ]: 20675 : for (auto p : m_symbcnodeset) {
617 [ + + ]: 61488 : for (const auto& s : g_cfg.get< tag::bc_sym >()) {
618 : : auto m = m_bnorm.find(s);
619 [ + + ]: 41248 : if (m != end(m_bnorm)) {
620 : : auto r = m->second.find(p);
621 [ + + ]: 36127 : if (r != end(m->second)) {
622 [ + - ]: 21515 : m_symbcnodes.push_back( p );
623 [ + - ]: 21515 : m_symbcnorms.push_back( r->second[0] );
624 [ + - ]: 21515 : m_symbcnorms.push_back( r->second[1] );
625 [ + - ]: 21515 : 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 [ + + ]: 1159219 : for (const auto& [p,q] : tk::lpoed) {
691 [ + + ]: 1105500 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
692 [ + + ]: 1630327 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
693 : 1105500 : d[f] = m_domedgeint.find( ed );
694 [ + + ]: 1105500 : if (d[f] == end(m_domedgeint)) break; else ++f;
695 : : }
696 [ + + ]: 392691 : if (f == 6) {
697 [ + - ]: 53719 : m_dsupedge[0].push_back( N[0] );
698 [ + - ]: 53719 : m_dsupedge[0].push_back( N[1] );
699 [ + - ]: 53719 : m_dsupedge[0].push_back( N[2] );
700 [ + - ]: 53719 : m_dsupedge[0].push_back( N[3] );
701 [ + + ]: 268595 : for (const auto& [a,b,c] : tk::lpofa) untri.erase( { N[a], N[b], N[c] } );
702 [ + + ]: 376033 : for (int ed=0; ed<6; ++ed) {
703 : : const auto& ded = d[ed]->second;
704 [ + - ]: 322314 : m_dsupint[0].push_back( sig[ed] * ded[0] );
705 [ + - ]: 322314 : m_dsupint[0].push_back( sig[ed] * ded[1] );
706 [ + - ][ + - ]: 322314 : m_dsupint[0].push_back( sig[ed] * ded[2] );
707 [ + - ]: 322314 : m_dsupint[0].push_back( ded[3] );
708 : 322314 : m_domedgeint.erase( d[ed] );
709 : : }
710 : : }
711 : : }
712 : :
713 [ + + ]: 623191 : for (const auto& N : untri) {
714 : : int f = 0;
715 : : tk::real sig[3];
716 [ + + ]: 2491024 : decltype(m_domedgeint)::const_iterator d[3];
717 [ + + ]: 987731 : for (const auto& [p,q] : tk::lpoet) {
718 [ + + ]: 957123 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
719 [ + + ]: 1464276 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
720 : 957123 : d[f] = m_domedgeint.find( ed );
721 [ + + ]: 957123 : if (d[f] == end(m_domedgeint)) break; else ++f;
722 : : }
723 [ + + ]: 622756 : if (f == 3) {
724 [ + - ]: 30608 : m_dsupedge[1].push_back( N[0] );
725 [ + - ]: 30608 : m_dsupedge[1].push_back( N[1] );
726 [ + - ]: 30608 : m_dsupedge[1].push_back( N[2] );
727 [ + + ]: 122432 : for (int ed=0; ed<3; ++ed) {
728 : : const auto& ded = d[ed]->second;
729 [ + - ]: 91824 : m_dsupint[1].push_back( sig[ed] * ded[0] );
730 [ + - ]: 91824 : m_dsupint[1].push_back( sig[ed] * ded[1] );
731 [ + - ][ + - ]: 91824 : m_dsupint[1].push_back( sig[ed] * ded[2] );
732 [ + - ]: 91824 : m_dsupint[1].push_back( ded[3] );
733 : 91824 : 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 [ + + ]: 127362 : for (const auto& [ed,d] : m_domedgeint) {
742 : 126927 : auto e = m_dsupedge[2].data() + k*2;
743 : 126927 : e[0] = tk::cref_find( lid, ed[0] );
744 : 126927 : e[1] = tk::cref_find( lid, ed[1] );
745 : 126927 : auto n = m_dsupint[2].data() + k*4;
746 : 126927 : n[0] = d[0];
747 : 126927 : n[1] = d[1];
748 : 126927 : n[2] = d[2];
749 : 126927 : n[3] = d[3];
750 : 126927 : ++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 [ + + ]: 8508 : for (auto g : nodes) {
805 : 8372 : auto i = tk::cref_find(lid,g);
806 [ + + ]: 89358 : for (auto e : tk::Around(esup,i)) {
807 : 80986 : const auto N = inpoel.data() + e*4;
808 [ + + ]: 566902 : for (const auto& [p,q] : tk::lpoed) {
809 : 485916 : tk::UnsMesh::Edge ged{ gid[N[p]], gid[N[q]] };
810 [ + - ]: 485916 : edges[ { N[p], N[q] } ] = tk::cref_find(m_domedgeint,ged)[3];
811 : : }
812 : : }
813 : : }
814 : : auto& ed = m_chbndedge[c];
815 [ + - ][ + + ]: 72203 : for (const auto& [e,sint] : edges) ed.emplace_back( e, sint );
816 : : }
817 : 19 : }
818 : :
819 : : void
820 : 5765 : ZalCG::deavol()
821 : : // *****************************************************************************
822 : : // Adjust node volumes along inactive neighbor chares
823 : : // *****************************************************************************
824 : : {
825 : 5765 : auto d = Disc();
826 [ + + ][ + - ]: 5765 : if (not g_cfg.get< tag::deactivate >() or d->NodeCommMap().empty()) return;
827 : :
828 : 259 : m_vol = d->Vol();
829 : : const auto& cvolc = d->Cvolc();
830 [ + + ]: 944 : for (auto i : m_inactive) {
831 [ + + ]: 33735 : for (const auto& [p,v] : tk::cref_find(cvolc,i)) {
832 : 33050 : 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 : 5765 : 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 : 5765 : auto d = Disc();
872 : :
873 : : // Apply Dirichlet BCs
874 [ + - ]: 5765 : physics::dirbc( u, t, d->Coord(), d->BoxNodes(), m_dirbcmasks );
875 : :
876 : : // Apply symmetry BCs
877 : 5765 : physics::symbc( u, m_symbcnodes, m_symbcnorms, /*pos=*/1 );
878 : :
879 : : // Apply farfield BCs
880 : 5765 : physics::farbc( u, m_farbcnodes, m_farbcnorms );
881 : :
882 : : // Apply pressure BCs
883 : 5765 : physics::prebc( u, m_prebcnodes, m_prebcvals );
884 : 5765 : }
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 : 5330 : deavol();
902 : :
903 : : // use constant dt if configured
904 [ + + ]: 5330 : 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 : 2410 : auto cfl = g_cfg.get< tag::cfl >();
913 : :
914 [ + + ]: 2410 : 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 [ + + ]: 832550 : for (std::size_t p=0; p<m_u.nunk(); ++p) {
932 : 830620 : auto r = m_u(p,0);
933 : 830620 : auto u = m_u(p,1)/r;
934 : 830620 : auto v = m_u(p,2)/r;
935 : 830620 : auto w = m_u(p,3)/r;
936 [ - + ]: 830620 : auto pr = eos::pressure( m_u(p,4) - 0.5*r*(u*u + v*v + w*w) );
937 [ - + ][ - + ]: 830620 : auto c = eos::soundspeed( r, std::max(pr,0.0) );
938 : 830620 : auto L = std::cbrt( vol[p] );
939 : 830620 : auto vel = std::sqrt( u*u + v*v + w*w );
940 [ - + ][ + + ]: 830620 : auto euler_dt = L / std::max( vel+c, 1.0e-8 );
941 : 830620 : mindt = std::min( mindt, euler_dt );
942 : : }
943 : 1930 : mindt *= cfl;
944 : :
945 : : }
946 : :
947 : 2410 : 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 : : // Detect blowup
975 : : auto eps = std::numeric_limits< tk::real >::epsilon();
976 [ - + ]: 6040 : if (newdt < eps) m_finished = 1;
977 : :
978 : : // Set new time step size
979 : 6040 : Disc()->setdt( newdt );
980 : :
981 [ + + ]: 6040 : if (m_deactivated) solve(); else rhs();
982 : 6040 : }
983 : :
984 : : void
985 : 5330 : ZalCG::rhs()
986 : : // *****************************************************************************
987 : : // Compute right-hand side of transport equations
988 : : // *****************************************************************************
989 : : {
990 : 5330 : auto d = Disc();
991 : : const auto& lid = d->Lid();
992 : :
993 : : // Compute own portion of right-hand side for all equations
994 : 5330 : zalesak::rhs( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_besym,
995 : 5330 : d->T(), d->Dt(), m_tp, m_dtp, m_u, m_rhs );
996 : :
997 : : // Communicate rhs to other chares on chare-boundary
998 [ + + ][ + + ]: 5330 : if (d->NodeCommMap().empty() or m_inactive.size() == d->NodeCommMap().size()){
999 : 145 : comrhs_complete();
1000 : : } else {
1001 [ + + ]: 51910 : for (const auto& [c,n] : d->NodeCommMap()) {
1002 : 665 : if (m_inactive.count(c)) continue;
1003 : : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
1004 [ + - ][ + + ]: 629390 : for (auto g : n) exp[g] = m_rhs[ tk::cref_find(lid,g) ];
1005 [ + - ][ + - ]: 92120 : thisProxy[c].comrhs( exp );
1006 : : }
1007 : : }
1008 : 5330 : ownrhs_complete();
1009 : 5330 : }
1010 : :
1011 : : void
1012 : 46060 : ZalCG::comrhs(
1013 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& inrhs )
1014 : : // *****************************************************************************
1015 : : // Receive contributions to right-hand side vector on chare-boundaries
1016 : : //! \param[in] inrhs Partial contributions of RHS to chare-boundary nodes. Key:
1017 : : //! global mesh node IDs, value: contributions for all scalar components.
1018 : : // *****************************************************************************
1019 : : {
1020 : : using tk::operator+=;
1021 [ + + ]: 1212720 : for (const auto& [g,r] : inrhs) m_rhsc[g] += r;
1022 : :
1023 : : // When we have heard from all chares we communicate with, this chare is done
1024 [ + + ]: 46060 : if (++m_nrhs + m_inactive.size() == Disc()->NodeCommMap().size()) {
1025 : 5185 : m_nrhs = 0;
1026 : 5185 : comrhs_complete();
1027 : : }
1028 : 46060 : }
1029 : :
1030 : : void
1031 : 5330 : ZalCG::fct()
1032 : : // *****************************************************************************
1033 : : // Continue with flux-corrected transport if enabled
1034 : : // *****************************************************************************
1035 : : {
1036 : 5330 : auto d = Disc();
1037 : : const auto& lid = d->Lid();
1038 : :
1039 : : // Combine own and communicated contributions to rhs
1040 [ + + ]: 442320 : for (const auto& [g,r] : m_rhsc) {
1041 : 436990 : auto i = tk::cref_find( lid, g );
1042 [ + + ]: 2657080 : for (std::size_t c=0; c<r.size(); ++c) m_rhs(i,c) += r[c];
1043 : : }
1044 : 5330 : tk::destroy(m_rhsc);
1045 : :
1046 [ + - ]: 5330 : if (g_cfg.get< tag::fct >()) aec(); else solve();
1047 : 5330 : }
1048 : :
1049 : : void
1050 : : // cppcheck-suppress unusedFunction
1051 : 5330 : ZalCG::aec()
1052 : : // *****************************************************************************
1053 : : // Compute antidiffusive contributions: P+/-
1054 : : // *****************************************************************************
1055 : : {
1056 : 5330 : auto d = Disc();
1057 : : const auto ncomp = m_u.nprop();
1058 : : const auto& lid = d->Lid();
1059 : :
1060 : : // Antidiffusive contributions: P+/-
1061 : :
1062 : 5330 : auto ctau = g_cfg.get< tag::fctdif >();
1063 : : m_p.fill( 0.0 );
1064 : :
1065 : : // tetrahedron superedges
1066 [ + + ]: 983015 : for (std::size_t e=0; e<m_dsupedge[0].size()/4; ++e) {
1067 : 977685 : const auto N = m_dsupedge[0].data() + e*4;
1068 : : const auto D = m_dsupint[0].data();
1069 : : std::size_t i = 0;
1070 [ + + ]: 6843795 : for (const auto& [p,q] : tk::lpoed) {
1071 : 5866110 : auto dif = D[(e*6+i)*4+3];
1072 [ + + ]: 35393940 : for (std::size_t c=0; c<ncomp; ++c) {
1073 [ + + ]: 29527830 : auto aec = -dif * ctau * (m_u(N[p],c) - m_u(N[q],c));
1074 : 29527830 : auto a = c*2;
1075 : 29527830 : auto b = a+1;
1076 [ + + ]: 29527830 : if (aec > 0.0) std::swap(a,b);
1077 : 29527830 : m_p(N[p],a) -= aec;
1078 : 29527830 : m_p(N[q],b) += aec;
1079 : : }
1080 : 5866110 : ++i;
1081 : : }
1082 : : }
1083 : :
1084 : : // triangle superedges
1085 [ + + ]: 536260 : for (std::size_t e=0; e<m_dsupedge[1].size()/3; ++e) {
1086 : 530930 : const auto N = m_dsupedge[1].data() + e*3;
1087 : : const auto D = m_dsupint[1].data();
1088 : : std::size_t i = 0;
1089 [ + + ]: 2123720 : for (const auto& [p,q] : tk::lpoet) {
1090 : 1592790 : auto dif = D[(e*3+i)*4+3];
1091 [ + + ]: 9650820 : for (std::size_t c=0; c<ncomp; ++c) {
1092 [ + + ]: 8058030 : auto aec = -dif * ctau * (m_u(N[p],c) - m_u(N[q],c));
1093 : 8058030 : auto a = c*2;
1094 : 8058030 : auto b = a+1;
1095 [ + + ]: 8058030 : if (aec > 0.0) std::swap(a,b);
1096 : 8058030 : m_p(N[p],a) -= aec;
1097 : 8058030 : m_p(N[q],b) += aec;
1098 : : }
1099 : 1592790 : ++i;
1100 : : }
1101 : : }
1102 : :
1103 : : // edges
1104 [ + + ]: 2276010 : for (std::size_t e=0; e<m_dsupedge[2].size()/2; ++e) {
1105 : 2270680 : const auto N = m_dsupedge[2].data() + e*2;
1106 : 2270680 : const auto dif = m_dsupint[2][e*4+3];
1107 [ + + ]: 13730700 : for (std::size_t c=0; c<ncomp; ++c) {
1108 [ + + ]: 11460020 : auto aec = -dif * ctau * (m_u(N[0],c) - m_u(N[1],c));
1109 : 11460020 : auto a = c*2;
1110 : 11460020 : auto b = a+1;
1111 [ + + ]: 11460020 : if (aec > 0.0) std::swap(a,b);
1112 : 11460020 : m_p(N[0],a) -= aec;
1113 : 11460020 : m_p(N[1],b) += aec;
1114 : : }
1115 : : }
1116 : :
1117 : : // Apply symmetry BCs on AEC
1118 [ + + ]: 376570 : for (std::size_t i=0; i<m_symbcnodes.size(); ++i) {
1119 : 371240 : auto p = m_symbcnodes[i];
1120 : 371240 : auto nx = m_symbcnorms[i*3+0];
1121 : 371240 : auto ny = m_symbcnorms[i*3+1];
1122 : 371240 : auto nz = m_symbcnorms[i*3+2];
1123 : 371240 : auto rvnp = m_p(p,2)*nx + m_p(p,4)*ny + m_p(p,6)*nz;
1124 : 371240 : auto rvnn = m_p(p,3)*nx + m_p(p,5)*ny + m_p(p,7)*nz;
1125 : 371240 : m_p(p,2) -= rvnp * nx;
1126 : 371240 : m_p(p,3) -= rvnn * nx;
1127 : 371240 : m_p(p,4) -= rvnp * ny;
1128 : 371240 : m_p(p,5) -= rvnn * ny;
1129 : 371240 : m_p(p,6) -= rvnp * nz;
1130 : 371240 : m_p(p,7) -= rvnn * nz;
1131 : : }
1132 : :
1133 : : // Communicate antidiffusive edge and low-order solution contributions
1134 [ + + ][ + + ]: 5330 : if (d->NodeCommMap().empty() or m_inactive.size() == d->NodeCommMap().size()){
1135 : 145 : comaec_complete();
1136 : : } else {
1137 [ + + ]: 51910 : for (const auto& [c,n] : d->NodeCommMap()) {
1138 : 665 : if (m_inactive.count(c)) continue;
1139 : : decltype(m_pc) exp;
1140 [ + - ][ + + ]: 629390 : for (auto g : n) exp[g] = m_p[ tk::cref_find(lid,g) ];
1141 [ + - ][ + - ]: 92120 : thisProxy[c].comaec( exp );
1142 : : }
1143 : : }
1144 : 5330 : ownaec_complete();
1145 : 5330 : }
1146 : :
1147 : : void
1148 : 46060 : ZalCG::comaec( const std::unordered_map< std::size_t,
1149 : : std::vector< tk::real > >& inaec )
1150 : : // *****************************************************************************
1151 : : // Receive antidiffusive and low-order contributions on chare-boundaries
1152 : : //! \param[in] inaec Partial contributions of antidiffusive edge and low-order
1153 : : //! solution contributions on chare-boundary nodes. Key: global mesh node IDs,
1154 : : //! value: 0: antidiffusive contributions, 1: low-order solution.
1155 : : // *****************************************************************************
1156 : : {
1157 : : using tk::operator+=;
1158 [ + + ]: 1212720 : for (const auto& [g,a] : inaec) m_pc[g] += a;
1159 : :
1160 : : // When we have heard from all chares we communicate with, this chare is done
1161 [ + + ]: 46060 : if (++m_naec + m_inactive.size() == Disc()->NodeCommMap().size()) {
1162 : 5185 : m_naec = 0;
1163 : 5185 : comaec_complete();
1164 : : }
1165 : 46060 : }
1166 : :
1167 : : void
1168 : 5330 : ZalCG::alw()
1169 : : // *****************************************************************************
1170 : : // Compute allowed limits, Q+/-
1171 : : // *****************************************************************************
1172 : : {
1173 : 5330 : auto d = Disc();
1174 : 5330 : const auto steady = g_cfg.get< tag::steady >();
1175 : : const auto npoin = m_u.nunk();
1176 : : const auto ncomp = m_u.nprop();
1177 : : const auto& lid = d->Lid();
1178 : :
1179 : : // Combine own and communicated contributions to antidiffusive contributions
1180 : : // and low-order solution
1181 [ + + ]: 442320 : for (const auto& [g,p] : m_pc) {
1182 : 436990 : auto i = tk::cref_find( lid, g );
1183 [ + + ]: 4877170 : for (std::size_t c=0; c<p.size(); ++c) m_p(i,c) += p[c];
1184 : : }
1185 : 5330 : tk::destroy(m_pc);
1186 : :
1187 : : // Finish computing antidiffusive contributions and low-order solution
1188 : : auto dt = d->Dt();
1189 [ + + ]: 1705230 : for (std::size_t i=0; i<npoin; ++i) {
1190 [ + + ]: 1699900 : if (steady) dt = m_dtp[i];
1191 [ + + ]: 10290420 : for (std::size_t c=0; c<ncomp; ++c) {
1192 : 8590520 : auto a = c*2;
1193 : 8590520 : auto b = a+1;
1194 : 8590520 : m_p(i,a) /= m_vol[i];
1195 : 8590520 : m_p(i,b) /= m_vol[i];
1196 : : // low-order solution
1197 : 8590520 : m_rhs(i,c) = m_u(i,c) - dt*m_rhs(i,c)/m_vol[i] - m_p(i,a) - m_p(i,b);
1198 : : }
1199 : : }
1200 : :
1201 : : // Allowed limits: Q+/-
1202 : :
1203 : : using std::max;
1204 : : using std::min;
1205 : :
1206 : : auto large = std::numeric_limits< tk::real >::max();
1207 [ + + ]: 1705230 : for (std::size_t i=0; i<m_q.nunk(); ++i) {
1208 [ + + ]: 10290420 : for (std::size_t c=0; c<m_q.nprop()/2; ++c) {
1209 : 8590520 : m_q(i,c*2+0) = -large;
1210 : 8590520 : m_q(i,c*2+1) = +large;
1211 : : }
1212 : : }
1213 : :
1214 : : // tetrahedron superedges
1215 [ + + ]: 983015 : for (std::size_t e=0; e<m_dsupedge[0].size()/4; ++e) {
1216 : 977685 : const auto N = m_dsupedge[0].data() + e*4;
1217 [ + + ]: 5898990 : for (std::size_t c=0; c<ncomp; ++c) {
1218 : 4921305 : auto a = c*2;
1219 : 4921305 : auto b = a+1;
1220 [ + + ]: 34449135 : for (const auto& [p,q] : tk::lpoed) {
1221 : : tk::real alwp, alwn;
1222 [ + + ]: 29527830 : if (g_cfg.get< tag::fctclip >()) {
1223 [ + + ][ + + ]: 648301 : alwp = max( m_rhs(N[p],c), m_rhs(N[q],c) );
1224 : 472500 : alwn = min( m_rhs(N[p],c), m_rhs(N[q],c) );
1225 : : } else {
1226 [ + + ][ + + ]: 58110660 : alwp = max( max(m_rhs(N[p],c), m_u(N[p],c)),
1227 [ + + ]: 29055330 : max(m_rhs(N[q],c), m_u(N[q],c)) );
1228 : 29055330 : alwn = min( min(m_rhs(N[p],c), m_u(N[p],c)),
1229 : : min(m_rhs(N[q],c), m_u(N[q],c)) );
1230 : : }
1231 [ + + ][ + + ]: 35761750 : m_q(N[p],a) = max(m_q(N[p],a), alwp);
1232 : 29527830 : m_q(N[p],b) = min(m_q(N[p],b), alwn);
1233 [ + + ][ + + ]: 39685757 : m_q(N[q],a) = max(m_q(N[q],a), alwp);
1234 : 29527830 : m_q(N[q],b) = min(m_q(N[q],b), alwn);
1235 : : }
1236 : : }
1237 : : }
1238 : :
1239 : : // triangle superedges
1240 [ + + ]: 536260 : for (std::size_t e=0; e<m_dsupedge[1].size()/3; ++e) {
1241 : 530930 : const auto N = m_dsupedge[1].data() + e*3;
1242 [ + + ]: 3216940 : for (std::size_t c=0; c<ncomp; ++c) {
1243 : 2686010 : auto a = c*2;
1244 : 2686010 : auto b = a+1;
1245 [ + + ]: 10744040 : for (const auto& [p,q] : tk::lpoet) {
1246 : : tk::real alwp, alwn;
1247 [ + + ]: 8058030 : if (g_cfg.get< tag::fctclip >()) {
1248 [ + + ][ + + ]: 254991 : alwp = max( m_rhs(N[p],c), m_rhs(N[q],c) );
1249 : 185100 : alwn = min( m_rhs(N[p],c), m_rhs(N[q],c) );
1250 : : } else {
1251 [ + + ][ + + ]: 15745860 : alwp = max( max(m_rhs(N[p],c), m_u(N[p],c)),
1252 [ + + ]: 7872930 : max(m_rhs(N[q],c), m_u(N[q],c)) );
1253 : 7872930 : alwn = min( min(m_rhs(N[p],c), m_u(N[p],c)),
1254 : : min(m_rhs(N[q],c), m_u(N[q],c)) );
1255 : : }
1256 [ + + ][ + + ]: 8996275 : m_q(N[p],a) = max(m_q(N[p],a), alwp);
1257 : 8058030 : m_q(N[p],b) = min(m_q(N[p],b), alwn);
1258 [ + + ][ + + ]: 9136151 : m_q(N[q],a) = max(m_q(N[q],a), alwp);
1259 : 8058030 : m_q(N[q],b) = min(m_q(N[q],b), alwn);
1260 : : }
1261 : : }
1262 : : }
1263 : :
1264 : : // edges
1265 [ + + ]: 2276010 : for (std::size_t e=0; e<m_dsupedge[2].size()/2; ++e) {
1266 : 2270680 : const auto N = m_dsupedge[2].data() + e*2;
1267 [ + + ]: 13730700 : for (std::size_t c=0; c<ncomp; ++c) {
1268 : 11460020 : auto a = c*2;
1269 : 11460020 : auto b = a+1;
1270 : : tk::real alwp, alwn;
1271 [ + + ]: 11460020 : if (g_cfg.get< tag::fctclip >()) {
1272 [ + + ][ + + ]: 294388 : alwp = max( m_rhs(N[0],c), m_rhs(N[1],c) );
1273 : 214000 : alwn = min( m_rhs(N[0],c), m_rhs(N[1],c) );
1274 : : } else {
1275 [ + + ][ + + ]: 22492040 : alwp = max( max(m_rhs(N[0],c), m_u(N[0],c)),
1276 [ + + ]: 11246020 : max(m_rhs(N[1],c), m_u(N[1],c)) );
1277 : 11246020 : alwn = min( min(m_rhs(N[0],c), m_u(N[0],c)),
1278 : : min(m_rhs(N[1],c), m_u(N[1],c)) );
1279 : : }
1280 [ + + ][ + + ]: 12427392 : m_q(N[0],a) = max(m_q(N[0],a), alwp);
1281 : 11460020 : m_q(N[0],b) = min(m_q(N[0],b), alwn);
1282 [ + + ][ + + ]: 12332793 : m_q(N[1],a) = max(m_q(N[1],a), alwp);
1283 : 11460020 : m_q(N[1],b) = min(m_q(N[1],b), alwn);
1284 : : }
1285 : : }
1286 : :
1287 : : // Communicate allowed limits contributions
1288 [ + + ][ + + ]: 5330 : if (d->NodeCommMap().empty() or m_inactive.size() == d->NodeCommMap().size()){
1289 : 145 : comalw_complete();
1290 : : } else {
1291 [ + + ]: 51910 : for (const auto& [c,n] : d->NodeCommMap()) {
1292 : 665 : if (m_inactive.count(c)) continue;
1293 : : decltype(m_qc) exp;
1294 [ + - ][ + + ]: 629390 : for (auto g : n) exp[g] = m_q[ tk::cref_find(lid,g) ];
1295 [ + - ][ + - ]: 92120 : thisProxy[c].comalw( exp );
1296 : : }
1297 : : }
1298 : 5330 : ownalw_complete();
1299 : 5330 : }
1300 : :
1301 : : void
1302 : 46060 : ZalCG::comalw( const std::unordered_map< std::size_t,
1303 : : std::vector< tk::real > >& inalw )
1304 : : // *****************************************************************************
1305 : : // Receive allowed limits contributions on chare-boundaries
1306 : : //! \param[in] inalw Partial contributions of allowed limits contributions on
1307 : : //! chare-boundary nodes. Key: global mesh node IDs, value: allowed limit
1308 : : //! contributions.
1309 : : // *****************************************************************************
1310 : : {
1311 [ + + ]: 629390 : for (const auto& [g,alw] : inalw) {
1312 : : auto& q = m_qc[g];
1313 : 583330 : q.resize( alw.size() );
1314 [ + + ]: 3544380 : for (std::size_t c=0; c<alw.size()/2; ++c) {
1315 : 2961050 : auto a = c*2;
1316 [ + + ]: 2961050 : auto b = a+1;
1317 [ + + ]: 2961050 : q[a] = std::max( q[a], alw[a] );
1318 : 2961050 : q[b] = std::min( q[b], alw[b] );
1319 : : }
1320 : : }
1321 : :
1322 : : // When we have heard from all chares we communicate with, this chare is done
1323 [ + + ]: 46060 : if (++m_nalw + m_inactive.size() == Disc()->NodeCommMap().size()) {
1324 : 5185 : m_nalw = 0;
1325 : 5185 : comalw_complete();
1326 : : }
1327 : 46060 : }
1328 : :
1329 : : void
1330 : 5330 : ZalCG::lim()
1331 : : // *****************************************************************************
1332 : : // Compute limit coefficients
1333 : : // *****************************************************************************
1334 : : {
1335 : 5330 : auto d = Disc();
1336 : : const auto npoin = m_u.nunk();
1337 : : const auto ncomp = m_u.nprop();
1338 : : const auto& lid = d->Lid();
1339 : :
1340 : : using std::max;
1341 : : using std::min;
1342 : :
1343 : : // Combine own and communicated contributions to allowed limits
1344 [ + + ]: 442320 : for (const auto& [g,alw] : m_qc) {
1345 : 436990 : auto i = tk::cref_find( lid, g );
1346 [ + + ]: 2657080 : for (std::size_t c=0; c<alw.size()/2; ++c) {
1347 : 2220090 : auto a = c*2;
1348 [ + + ]: 2220090 : auto b = a+1;
1349 [ + + ]: 2220090 : m_q(i,a) = max( m_q(i,a), alw[a] );
1350 : 2220090 : m_q(i,b) = min( m_q(i,b), alw[b] );
1351 : : }
1352 : : }
1353 : 5330 : tk::destroy(m_qc);
1354 : :
1355 : : // Finish computing allowed limits
1356 [ + + ]: 1705230 : for (std::size_t i=0; i<npoin; ++i) {
1357 [ + + ]: 10290420 : for (std::size_t c=0; c<ncomp; ++c) {
1358 : 8590520 : auto a = c*2;
1359 : 8590520 : auto b = a+1;
1360 : 8590520 : m_q(i,a) -= m_rhs(i,c);
1361 : 8590520 : m_q(i,b) -= m_rhs(i,c);
1362 : : }
1363 : : }
1364 : :
1365 : : // Limit coefficients, C
1366 : :
1367 [ + + ]: 1705230 : for (std::size_t i=0; i<npoin; ++i) {
1368 [ + + ]: 10290420 : for (std::size_t c=0; c<ncomp; ++c) {
1369 : 8590520 : auto a = c*2;
1370 [ + + ]: 8590520 : auto b = a+1;
1371 : : auto eps = std::numeric_limits< tk::real >::epsilon();
1372 [ + + ][ + + ]: 8643661 : m_q(i,a) = m_p(i,a) < eps ? 0.0 : min(1.0, m_q(i,a)/m_p(i,a));
1373 [ + + ][ + + ]: 8644097 : m_q(i,b) = m_p(i,b) > -eps ? 0.0 : min(1.0, m_q(i,b)/m_p(i,b));
1374 : : }
1375 : : }
1376 : :
1377 : : // Limited antidiffusive contributions
1378 : :
1379 : 5330 : auto ctau = g_cfg.get< tag::fctdif >();
1380 : : m_a.fill( 0.0 );
1381 : :
1382 : 5330 : auto fctsys = g_cfg.get< tag::fctsys >();
1383 [ + + ]: 9310 : for (auto& c : fctsys) --c;
1384 : :
1385 : : #if defined(__clang__)
1386 : : #pragma clang diagnostic push
1387 : : #pragma clang diagnostic ignored "-Wvla"
1388 : : #pragma clang diagnostic ignored "-Wvla-extension"
1389 : : #elif defined(STRICT_GNUC)
1390 : : #pragma GCC diagnostic push
1391 : : #pragma GCC diagnostic ignored "-Wvla"
1392 : : #endif
1393 : :
1394 : : // tetrahedron superedges
1395 [ + + ]: 983015 : for (std::size_t e=0; e<m_dsupedge[0].size()/4; ++e) {
1396 : 977685 : const auto N = m_dsupedge[0].data() + e*4;
1397 : : const auto D = m_dsupint[0].data();
1398 : : auto C = m_dsuplim[0].data();
1399 : : std::size_t i = 0;
1400 [ + + ]: 6843795 : for (const auto& [p,q] : tk::lpoed) {
1401 : 5866110 : auto dif = D[(e*6+i)*4+3];
1402 : 5866110 : auto coef = C + (e*6+i)*ncomp;
1403 : 5866110 : tk::real aec[ncomp];
1404 [ + + ]: 35393940 : for (std::size_t c=0; c<ncomp; ++c) {
1405 [ + + ]: 29527830 : aec[c] = -dif * ctau * (m_u(N[p],c) - m_u(N[q],c));
1406 [ + + ]: 29527830 : if (m_fctfreeze) continue;
1407 : 25777530 : auto a = c*2;
1408 : 25777530 : auto b = a+1;
1409 [ + + ][ + + ]: 77332590 : coef[c] = min( aec[c] < 0.0 ? m_q(N[p],a) : m_q(N[p],b),
1410 : : aec[c] > 0.0 ? m_q(N[q],a) : m_q(N[q],b) );
1411 : : }
1412 : 5866110 : tk::real cs = 1.0;
1413 [ + + ][ + + ]: 18033360 : for (auto c : fctsys) cs = min( cs, coef[c] );
1414 [ + + ]: 18033360 : for (auto c : fctsys) coef[c] = cs;
1415 [ + + ]: 35393940 : for (std::size_t c=0; c<ncomp; ++c) {
1416 : 29527830 : aec[c] *= coef[c];
1417 : 29527830 : m_a(N[p],c) -= aec[c];
1418 : 29527830 : m_a(N[q],c) += aec[c];
1419 : : }
1420 : 5866110 : ++i;
1421 : 5866110 : }
1422 : : }
1423 : :
1424 : : // triangle superedges
1425 [ + + ]: 536260 : for (std::size_t e=0; e<m_dsupedge[1].size()/3; ++e) {
1426 : 530930 : const auto N = m_dsupedge[1].data() + e*3;
1427 : : const auto D = m_dsupint[1].data();
1428 : : auto C = m_dsuplim[0].data();
1429 : : std::size_t i = 0;
1430 [ + + ]: 2123720 : for (const auto& [p,q] : tk::lpoet) {
1431 : 1592790 : auto dif = D[(e*3+i)*4+3];
1432 : 1592790 : auto coef = C + (e*3+i)*ncomp;
1433 : 1592790 : tk::real aec[ncomp];
1434 [ + + ]: 9650820 : for (std::size_t c=0; c<ncomp; ++c) {
1435 [ + + ]: 8058030 : aec[c] = -dif * ctau * (m_u(N[p],c) - m_u(N[q],c));
1436 [ + + ]: 8058030 : if (m_fctfreeze) continue;
1437 : 7300680 : auto a = c*2;
1438 : 7300680 : auto b = a+1;
1439 [ + + ][ + + ]: 21902040 : coef[c] = min( aec[c] < 0.0 ? m_q(N[p],a) : m_q(N[p],b),
1440 : : aec[c] > 0.0 ? m_q(N[q],a) : m_q(N[q],b) );
1441 : : }
1442 : 1592790 : tk::real cs = 1.0;
1443 [ + + ][ + + ]: 5521050 : for (auto c : fctsys) cs = min( cs, coef[c] );
1444 [ + + ]: 5521050 : for (auto c : fctsys) coef[c] = cs;
1445 [ + + ]: 9650820 : for (std::size_t c=0; c<ncomp; ++c) {
1446 : 8058030 : aec[c] *= coef[c];
1447 : 8058030 : m_a(N[p],c) -= aec[c];
1448 : 8058030 : m_a(N[q],c) += aec[c];
1449 : : }
1450 : 1592790 : ++i;
1451 : 1592790 : }
1452 : : }
1453 : :
1454 : : // edges
1455 [ + + ]: 2276010 : for (std::size_t e=0; e<m_dsupedge[2].size()/2; ++e) {
1456 : 2270680 : const auto N = m_dsupedge[2].data() + e*2;
1457 : 2270680 : const auto dif = m_dsupint[2][e*4+3];
1458 : 2270680 : auto coef = m_dsuplim[2].data() + e*ncomp;
1459 : 2270680 : tk::real aec[ncomp];
1460 [ + + ]: 13730700 : for (std::size_t c=0; c<ncomp; ++c) {
1461 [ + + ]: 11460020 : aec[c] = -dif * ctau * (m_u(N[0],c) - m_u(N[1],c));
1462 : 11460020 : auto a = c*2;
1463 : 11460020 : auto b = a+1;
1464 [ + + ]: 11460020 : if (m_fctfreeze) continue;
1465 [ + + ][ + + ]: 30574635 : coef[c] = min( aec[c] < 0.0 ? m_q(N[0],a) : m_q(N[0],b),
1466 : : aec[c] > 0.0 ? m_q(N[1],a) : m_q(N[1],b) );
1467 : : }
1468 : 2270680 : tk::real cs = 1.0;
1469 [ + + ][ + + ]: 7359530 : for (auto c : fctsys) cs = min( cs, coef[c] );
1470 [ + + ]: 7359530 : for (auto c : fctsys) coef[c] = cs;
1471 [ + + ]: 13730700 : for (std::size_t c=0; c<ncomp; ++c) {
1472 : 11460020 : aec[c] *= coef[c];
1473 : 11460020 : m_a(N[0],c) -= aec[c];
1474 : 11460020 : m_a(N[1],c) += aec[c];
1475 : : }
1476 : 2270680 : }
1477 : :
1478 : : #if defined(__clang__)
1479 : : #pragma clang diagnostic pop
1480 : : #elif defined(STRICT_GNUC)
1481 : : #pragma GCC diagnostic pop
1482 : : #endif
1483 : :
1484 : : // Communicate limited antidiffusive contributions
1485 [ + + ][ + + ]: 5330 : if (d->NodeCommMap().empty() or m_inactive.size() == d->NodeCommMap().size()){
1486 [ + - ]: 145 : comlim_complete();
1487 : : } else {
1488 [ + + ]: 51910 : for (const auto& [c,n] : d->NodeCommMap()) {
1489 : 665 : if (m_inactive.count(c)) continue;
1490 : : decltype(m_ac) exp;
1491 [ + - ][ + + ]: 629390 : for (auto g : n) exp[g] = m_a[ tk::cref_find(lid,g) ];
1492 [ + - ][ + - ]: 92120 : thisProxy[c].comlim( exp );
1493 : : }
1494 : : }
1495 [ + - ]: 5330 : ownlim_complete();
1496 : 5330 : }
1497 : :
1498 : : void
1499 : 46060 : ZalCG::comlim( const std::unordered_map< std::size_t,
1500 : : std::vector< tk::real > >& inlim )
1501 : : // *****************************************************************************
1502 : : // Receive limited antidiffusive contributions on chare-boundaries
1503 : : //! \param[in] inlim Partial contributions of limited contributions on
1504 : : //! chare-boundary nodes. Key: global mesh node IDs, value: limited
1505 : : //! contributions.
1506 : : // *****************************************************************************
1507 : : {
1508 : : using tk::operator+=;
1509 [ + + ]: 1212720 : for (const auto& [g,a] : inlim) m_ac[g] += a;
1510 : :
1511 : : // When we have heard from all chares we communicate with, this chare is done
1512 [ + + ]: 46060 : if (++m_nlim + m_inactive.size() == Disc()->NodeCommMap().size()) {
1513 : 5185 : m_nlim = 0;
1514 : 5185 : comlim_complete();
1515 : : }
1516 : 46060 : }
1517 : :
1518 : : void
1519 : 6040 : ZalCG::solve()
1520 : : // *****************************************************************************
1521 : : // Advance systems of equations
1522 : : // *****************************************************************************
1523 : : {
1524 : 6040 : auto d = Disc();
1525 : : const auto npoin = m_u.nunk();
1526 : : const auto ncomp = m_u.nprop();
1527 : 6040 : const auto steady = g_cfg.get< tag::steady >();
1528 : :
1529 [ + + ]: 6040 : if (m_deactivated) {
1530 : :
1531 : : m_a = m_u;
1532 : :
1533 : : } else {
1534 : :
1535 : : // Combine own and communicated contributions to limited antidiffusive
1536 : : // contributions
1537 [ + + ]: 442320 : for (const auto& [g,a] : m_ac) {
1538 : 436990 : auto i = tk::cref_find( d->Lid(), g );
1539 [ + + ]: 2657080 : for (std::size_t c=0; c<a.size(); ++c) m_a(i,c) += a[c];
1540 : : }
1541 : 5330 : tk::destroy(m_ac);
1542 : :
1543 : : tk::Fields u;
1544 [ + + ]: 5330 : std::size_t cstart = m_freezeflow > 1.0 ? 5 : 0;
1545 : : if (cstart) u = m_u;
1546 : :
1547 [ + - ]: 5330 : if (g_cfg.get< tag::fct >()) {
1548 : : // Apply limited antidiffusive contributions to low-order solution
1549 [ + + ]: 1705230 : for (std::size_t i=0; i<npoin; ++i) {
1550 [ + + ]: 10290420 : for (std::size_t c=0; c<ncomp; ++c) {
1551 : 8590520 : m_a(i,c) = m_rhs(i,c) + m_a(i,c)/m_vol[i];
1552 : : }
1553 : : }
1554 : : } else {
1555 : : // Apply rhs
1556 : : auto dt = d->Dt();
1557 [ - - ]: 0 : for (std::size_t i=0; i<npoin; ++i) {
1558 [ - - ]: 0 : if (steady) dt = m_dtp[i];
1559 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c) {
1560 : 0 : m_a(i,c) = m_u(i,c) - dt*m_rhs(i,c)/m_vol[i];
1561 : : }
1562 : : }
1563 : : }
1564 : :
1565 : : // Configure and apply scalar source to solution (if defined)
1566 [ + - ]: 5330 : auto src = problems::PHYS_SRC();
1567 [ - + ][ - - ]: 5330 : if (src) src( d->Coord(), d->T(), m_a );
1568 : :
1569 : : // Freeze flow if configured and apply multiplier on scalar(s)
1570 [ + + ]: 5330 : if (d->T() > g_cfg.get< tag::freezetime >()) {
1571 : 4895 : m_freezeflow = g_cfg.get< tag::freezeflow >();
1572 : : }
1573 : :
1574 : : // Enforce boundary conditions
1575 [ + - ]: 5330 : BC( m_a, d->T() + d->Dt() );
1576 : :
1577 : : // Explicitly zero out flow for freezeflow
1578 [ + + ]: 5330 : if (cstart) {
1579 [ + + ]: 82800 : for (std::size_t i=0; i<npoin; ++i) {
1580 [ + + ]: 491508 : for (std::size_t c=0; c<cstart; ++c) {
1581 : 409590 : m_a(i,c) = u(i,c);
1582 : : }
1583 : : }
1584 : : }
1585 : :
1586 : : }
1587 : :
1588 : : // Compute diagnostics, e.g., residuals
1589 : 6040 : auto diag_iter = g_cfg.get< tag::diag_iter >();
1590 : 6040 : auto diag = m_diag.rhocompute( *d, m_a, m_u, diag_iter );
1591 : :
1592 : : // Update solution
1593 : : m_u = m_a;
1594 : : m_a.fill( 0.0 );
1595 : :
1596 : : // Increase number of iterations and physical time
1597 : 6040 : d->next();
1598 : :
1599 : : // Advance physical time for local time stepping
1600 [ + + ]: 6040 : if (steady) {
1601 : : using tk::operator+=;
1602 : 480 : m_tp += m_dtp;
1603 : : }
1604 : :
1605 : : // Evaluate residuals
1606 [ - + ][ - - ]: 6040 : if (!diag) evalres( std::vector< tk::real >( ncomp, 1.0 ) );
1607 : 6040 : }
1608 : :
1609 : : void
1610 : 6040 : ZalCG::evalres( const std::vector< tk::real >& l2res )
1611 : : // *****************************************************************************
1612 : : // Evaluate residuals
1613 : : //! \param[in] l2res L2-norms of the residual for each scalar component
1614 : : //! computed across the whole problem
1615 : : // *****************************************************************************
1616 : : {
1617 : 6040 : auto d = Disc();
1618 : :
1619 [ + + ]: 6040 : if (g_cfg.get< tag::steady >()) {
1620 : 480 : const auto rc = g_cfg.get< tag::rescomp >() - 1;
1621 : 480 : d->residual( l2res[rc] );
1622 [ + + ]: 480 : if (l2res[rc] < g_cfg.get< tag::fctfreeze >()) m_fctfreeze = 1;
1623 : : }
1624 : :
1625 [ + + ]: 6040 : if (d->deactivate()) deactivate(); else refine();
1626 : 6040 : }
1627 : :
1628 : : int
1629 : 131307 : ZalCG::active( std::size_t p,
1630 : : std::size_t q,
1631 : : const std::vector< uint64_t >& sys ) const
1632 : : // *****************************************************************************
1633 : : // Decide if edge is active
1634 : : //! \param[in] p Local id of left edge-end point
1635 : : //! \param[in] q Local id of right edge-end point
1636 : : //! \param[in] sys List of components to consider
1637 : : //! \return True if active, false if not
1638 : : // *****************************************************************************
1639 : : {
1640 : 131307 : auto tol = g_cfg.get< tag::deatol >();
1641 : :
1642 : 131307 : tk::real maxdiff = 0.0;
1643 [ + + ][ + + ]: 262614 : for (auto c : sys) {
1644 [ + + ]: 131307 : auto e = std::abs( m_u(p,c) - m_u(q,c) );
1645 : 131307 : maxdiff = std::max( maxdiff, e );
1646 : : }
1647 : :
1648 : 131307 : return maxdiff > tol;
1649 : : }
1650 : :
1651 : : int
1652 : 19 : ZalCG::dea( const std::vector< uint64_t >& sys ) const
1653 : : // *****************************************************************************
1654 : : // Decide whether to deactivate this chare
1655 : : //! \param[in] sys List of components to consider
1656 : : //! \return Nonzero to deactivate, zero to keep active
1657 : : // *****************************************************************************
1658 : : {
1659 [ - + ]: 19 : if (m_toreactivate) return 0;
1660 [ - + ]: 19 : if (m_deactivated) return 1;
1661 : :
1662 : : // tetrahedron superedges
1663 [ + + ]: 7324 : for (std::size_t e=0; e<m_dsupedge[0].size()/4; ++e) {
1664 : 7306 : const auto N = m_dsupedge[0].data() + e*4;
1665 [ + + ]: 51137 : for (const auto& [p,q] : tk::lpoed) {
1666 [ + + ]: 43832 : if (active(N[p],N[q],sys)) return 0;
1667 : : }
1668 : : }
1669 : :
1670 : : // triangle superedges
1671 [ + + ]: 5270 : for (std::size_t e=0; e<m_dsupedge[1].size()/3; ++e) {
1672 : 5252 : const auto N = m_dsupedge[1].data() + e*3;
1673 [ + + ]: 21008 : for (const auto& [p,q] : tk::lpoet) {
1674 [ - + ]: 15756 : if (active(N[p],N[q],sys)) return 0;
1675 : : }
1676 : : }
1677 : :
1678 : : // edges
1679 [ + + ]: 19280 : for (std::size_t e=0; e<m_dsupedge[2].size()/2; ++e) {
1680 : 19262 : const auto N = m_dsupedge[2].data() + e*2;
1681 [ - + ]: 19262 : if (active(N[0],N[1],sys)) return 0;
1682 : : }
1683 : :
1684 : : return 1;
1685 : : }
1686 : :
1687 : : void
1688 : 48 : ZalCG::rea( const std::vector< uint64_t >& sys, std::unordered_set< int >& req )
1689 : : const
1690 : : // *****************************************************************************
1691 : : // Decide whether to reactivate any neighbor chare
1692 : : //! \param[in] sys List of components to consider
1693 : : //! \param[in,out] req Set of neighbor chares to reactivate
1694 : : // *****************************************************************************
1695 : : {
1696 [ + + ]: 401 : for (const auto& [c,edges] : m_chbndedge) {
1697 : 216 : if (not m_inactive.count(c)) continue;
1698 [ + + ]: 52588 : for (const auto& [e,sint] : edges) {
1699 [ + - ][ + + ]: 52457 : if (active(e[0],e[1],sys)) {
1700 : : req.insert(c);
1701 : : break;
1702 : : }
1703 : : }
1704 : : }
1705 : 48 : }
1706 : :
1707 : : void
1708 : 48 : ZalCG::huldif()
1709 : : // *****************************************************************************
1710 : : // Apply diffusion on active hull
1711 : : // *****************************************************************************
1712 : : {
1713 : : const auto eps = std::numeric_limits< tk::real >::epsilon();
1714 [ + - ]: 48 : const auto dif = g_cfg.get< tag::deadif >();
1715 [ + - ]: 48 : if (std::abs(dif) < eps) return;
1716 : :
1717 : : const auto npoin = m_u.nunk();
1718 : : const auto ncomp = m_u.nprop();
1719 : :
1720 : : m_rhs.fill( 0.0 );
1721 [ + + ]: 401 : for (const auto& [ch,edges] : m_chbndedge) {
1722 : 216 : if (not m_inactive.count(ch)) continue;
1723 [ + + ]: 56673 : for (const auto& [e,sint] : edges) {
1724 : 56536 : auto p = e[0];
1725 : 56536 : auto q = e[1];
1726 : 56536 : auto fw = dif * sint;
1727 [ + + ]: 339216 : for (std::size_t c=0; c<ncomp; ++c) {
1728 : 282680 : auto f = fw*(m_u(p,c) - m_u(q,c));
1729 : 282680 : m_rhs(p,c) -= f;
1730 : 282680 : m_rhs(q,c) += f;
1731 : : }
1732 : : }
1733 : : }
1734 : :
1735 [ + + ]: 37506 : for (std::size_t i=0; i<npoin; ++i) {
1736 [ + + ]: 224748 : for (std::size_t c=0; c<ncomp; ++c) {
1737 : 187290 : m_u(i,c) += m_rhs(i,c)/m_vol[i];
1738 : : }
1739 : : }
1740 : : }
1741 : :
1742 : : void
1743 : 190 : ZalCG::deactivate()
1744 : : // *****************************************************************************
1745 : : // Deactivate regions
1746 : : // *****************************************************************************
1747 : : {
1748 : 190 : auto d = Disc();
1749 : :
1750 : : std::unordered_set< int > reactivate;
1751 : :
1752 [ + + ]: 190 : if (not m_deactivated) {
1753 : : // Apply diffusion on active hull
1754 [ + - ]: 48 : huldif();
1755 : : // Query scalar components to deactivate based on
1756 [ + - ]: 48 : auto sys = g_cfg.get< tag::deasys >();
1757 [ + + ]: 96 : for (auto& c : sys) --c;
1758 : : // Decide whether to deactivate this chare
1759 [ + - ][ + + ]: 48 : if (d->deastart() and dea(sys)) m_todeactivate = 1;
[ + - ][ + + ]
1760 : : // Decide whether to reactivate any neighbor chare
1761 [ + - ]: 48 : rea( sys, reactivate );
1762 : : }
1763 : :
1764 : : // Communicate deactivation requests
1765 [ + + ]: 1550 : for (const auto& [c,n] : d->NodeCommMap()) {
1766 [ + - ][ + - ]: 2720 : thisProxy[c].comdea( reactivate.count(c) );
[ - + ][ - - ]
1767 : : }
1768 [ + - ]: 190 : owndea_complete();
1769 : 190 : }
1770 : :
1771 : : void
1772 : 1360 : ZalCG::comdea( std::size_t reactivate )
1773 : : // *****************************************************************************
1774 : : // Communicate deactivation desires
1775 : : //! \param[in] reactivate 1 if sender wants to reactivate this chare, 0 if not
1776 : : // *****************************************************************************
1777 : : {
1778 : 1360 : auto d = Disc();
1779 : :
1780 [ + + ][ + + ]: 1360 : if (m_deactivated and reactivate) m_toreactivate = 1;
1781 : :
1782 : : // Communicate activation status to neighbors
1783 [ + + ]: 1360 : if (++m_ndea == d->NodeCommMap().size()) {
1784 : 190 : m_ndea = 0;
1785 : 190 : comdea_complete();
1786 : : }
1787 : 1360 : }
1788 : :
1789 : : void
1790 : 190 : ZalCG::activate()
1791 : : // *****************************************************************************
1792 : : // Compute deactivation status
1793 : : // *****************************************************************************
1794 : : {
1795 : 190 : auto d = Disc();
1796 : :
1797 [ + + ]: 190 : if (m_todeactivate) m_deactivated = 1;
1798 [ + + ]: 190 : if (m_toreactivate) m_deactivated = 0;
1799 : :
1800 : : // Communicate deactivation status
1801 [ + + ]: 1550 : for (const auto& [c,n] : d->NodeCommMap()) {
1802 [ + - ]: 2720 : thisProxy[c].comact( thisIndex, m_deactivated );
1803 : : }
1804 : 190 : ownact_complete();
1805 : 190 : }
1806 : :
1807 : : void
1808 : 1360 : ZalCG::comact( int ch, int deactivated )
1809 : : // *****************************************************************************
1810 : : // Communicate deactivation status
1811 : : //! \param[in] ch Sender chare id
1812 : : //! \param[in] deactivated 1 if sender was deactivated, 0 if not
1813 : : // *****************************************************************************
1814 : : {
1815 : 1360 : auto d = Disc();
1816 : :
1817 [ + + ]: 1360 : if (deactivated) m_inactive.insert(ch); else m_inactive.erase(ch);
1818 : :
1819 : : // When we have heard from all chares we communicate with, this chare is done
1820 [ + + ]: 1360 : if (++m_nact == d->NodeCommMap().size()) {
1821 : : // Aggregate deactivation status to every chare
1822 [ + - ]: 190 : contribute( sizeof(int), &m_deactivated, CkReduction::sum_int,
1823 : 190 : CkCallback(CkReductionTarget(ZalCG,deastat), thisProxy) );
1824 : 190 : m_todeactivate = 0;
1825 : 190 : m_toreactivate = 0;
1826 : 190 : m_nact = 0;
1827 : 190 : comact_complete();
1828 : : }
1829 : 1360 : }
1830 : :
1831 : : void
1832 : 190 : ZalCG::deastat( int dea )
1833 : : // *****************************************************************************
1834 : : // Receive aggregate deactivation status
1835 : : //! \param[in] dea Sum of deactived chares
1836 : : // *****************************************************************************
1837 : : {
1838 : 190 : auto d = Disc();
1839 : :
1840 : : // Disallow deactivating all chares
1841 [ - + ]: 190 : if (dea == d->nchare()) {
1842 : : dea = 0;
1843 : 0 : m_deactivated = 0;
1844 : : m_inactive.clear();
1845 : : }
1846 : :
1847 : : // Report number of deactivated chares for diagnostics
1848 [ + + ]: 190 : if (thisIndex == 0) d->deactivated( dea );
1849 : :
1850 : 190 : deastat_complete();
1851 : 190 : }
1852 : :
1853 : : void
1854 : 6040 : ZalCG::refine()
1855 : : // *****************************************************************************
1856 : : // Optionally refine/derefine mesh
1857 : : // *****************************************************************************
1858 : : {
1859 : 6040 : auto d = Disc();
1860 : :
1861 : : // See if this is the last time step
1862 [ + + ]: 6040 : if (d->finished()) m_finished = 1;
1863 : :
1864 : 6040 : auto dtref = g_cfg.get< tag::href_dt >();
1865 : 6040 : auto dtfreq = g_cfg.get< tag::href_dtfreq >();
1866 : :
1867 : : // if t>0 refinement enabled and we hit the frequency
1868 [ - + ][ - - ]: 6040 : if (dtref && !(d->It() % dtfreq)) { // refine
1869 : :
1870 : 0 : d->refined() = 1;
1871 : 0 : d->startvol();
1872 : 0 : d->Ref()->dtref( m_bface, m_bnode, m_triinpoel );
1873 : :
1874 : : // Activate SDAG waits for re-computing the integrals
1875 [ - - ]: 0 : thisProxy[ thisIndex ].wait4int();
1876 : :
1877 : : } else { // do not refine
1878 : :
1879 : 6040 : d->refined() = 0;
1880 : 6040 : feop_complete();
1881 : 6040 : resize_complete();
1882 : :
1883 : : }
1884 : 6040 : }
1885 : :
1886 : : void
1887 : 0 : ZalCG::resizePostAMR(
1888 : : const std::vector< std::size_t >& /*ginpoel*/,
1889 : : const tk::UnsMesh::Chunk& chunk,
1890 : : const tk::UnsMesh::Coords& coord,
1891 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
1892 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
1893 : : const std::set< std::size_t >& removedNodes,
1894 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
1895 : : nodeCommMap,
1896 : : const std::map< int, std::vector< std::size_t > >& bface,
1897 : : const std::map< int, std::vector< std::size_t > >& bnode,
1898 : : const std::vector< std::size_t >& triinpoel )
1899 : : // *****************************************************************************
1900 : : // Receive new mesh from Refiner
1901 : : //! \param[in] ginpoel Mesh connectivity with global node ids
1902 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
1903 : : //! \param[in] coord New mesh node coordinates
1904 : : //! \param[in] addedNodes Newly added mesh nodes and their parents (local ids)
1905 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1906 : : //! \param[in] removedNodes Newly removed mesh node local ids
1907 : : //! \param[in] nodeCommMap New node communication map
1908 : : //! \param[in] bface Boundary-faces mapped to side set ids
1909 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
1910 : : //! \param[in] triinpoel Boundary-face connectivity
1911 : : // *****************************************************************************
1912 : : {
1913 : 0 : auto d = Disc();
1914 : : auto ncomp = m_u.nprop();
1915 : :
1916 : 0 : d->Itf() = 0; // Zero field output iteration count if AMR
1917 : 0 : ++d->Itr(); // Increase number of iterations with a change in the mesh
1918 : :
1919 : : // Resize mesh data structures after mesh refinement
1920 : 0 : d->resizePostAMR( chunk, coord, nodeCommMap, removedNodes );
1921 : :
1922 : : Assert(coord[0].size() == m_u.nunk()-removedNodes.size()+addedNodes.size(),
1923 : : "Incorrect vector length post-AMR: expected length after resizing = " +
1924 : : std::to_string(coord[0].size()) + ", actual unknown vector length = " +
1925 : : std::to_string(m_u.nunk()-removedNodes.size()+addedNodes.size()));
1926 : :
1927 : : // Remove newly removed nodes from solution vectors
1928 : 0 : m_u.rm( removedNodes );
1929 : 0 : m_rhs.rm( removedNodes );
1930 : :
1931 : : // Resize auxiliary solution vectors
1932 : : auto npoin = coord[0].size();
1933 : : m_u.resize( npoin );
1934 : : m_rhs.resize( npoin );
1935 : :
1936 : : // Update solution on new mesh
1937 [ - - ]: 0 : for (const auto& n : addedNodes)
1938 [ - - ]: 0 : for (std::size_t c=0; c<ncomp; ++c) {
1939 : : Assert(n.first < m_u.nunk(), "Added node index out of bounds post-AMR");
1940 : : Assert(n.second[0] < m_u.nunk() && n.second[1] < m_u.nunk(),
1941 : : "Indices of parent-edge nodes out of bounds post-AMR");
1942 : 0 : m_u(n.first,c) = (m_u(n.second[0],c) + m_u(n.second[1],c))/2.0;
1943 : : }
1944 : :
1945 : : // Update physical-boundary node-, face-, and element lists
1946 : : m_bnode = bnode;
1947 : : m_bface = bface;
1948 : 0 : m_triinpoel = tk::remap( triinpoel, d->Lid() );
1949 : :
1950 : 0 : auto meshid = d->MeshId();
1951 [ - - ]: 0 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1952 : 0 : CkCallback(CkReductionTarget(Transporter,resized), d->Tr()) );
1953 : 0 : }
1954 : :
1955 : : void
1956 : 1181 : ZalCG::writeFields( CkCallback cb )
1957 : : // *****************************************************************************
1958 : : // Output mesh-based fields to file
1959 : : //! \param[in] cb Function to continue with after the write
1960 : : // *****************************************************************************
1961 : : {
1962 [ + + ]: 1181 : if (g_cfg.get< tag::benchmark >()) { cb.send(); return; }
1963 : :
1964 : 597 : auto d = Disc();
1965 : : auto ncomp = m_u.nprop();
1966 : 597 : auto steady = g_cfg.get< tag::steady >();
1967 : :
1968 : : // Field output
1969 : :
1970 : : std::vector< std::string > nodefieldnames
1971 [ - + ][ + + ]: 4179 : {"density", "velocityx", "velocityy", "velocityz", "energy", "pressure"};
[ - + ][ - - ]
[ - - ]
1972 [ + - ]: 144 : if (steady) nodefieldnames.push_back( "mach" );
1973 : :
1974 : : using tk::operator/=;
1975 [ + - ]: 597 : auto r = m_u.extract(0);
1976 [ + - ][ + - ]: 597 : auto u = m_u.extract(1); u /= r;
1977 [ + - ][ + - ]: 597 : auto v = m_u.extract(2); v /= r;
1978 [ + - ][ + - ]: 597 : auto w = m_u.extract(3); w /= r;
1979 [ + - ][ + - ]: 597 : auto e = m_u.extract(4); e /= r;
1980 [ + - ][ + + ]: 597 : std::vector< tk::real > pr( m_u.nunk() ), ma;
[ - - ]
1981 [ + + ][ + - ]: 597 : if (steady) ma.resize( m_u.nunk() );
1982 [ + + ]: 437248 : for (std::size_t i=0; i<pr.size(); ++i) {
1983 : 436651 : auto vv = u[i]*u[i] + v[i]*v[i] + w[i]*w[i];
1984 [ + + ]: 436651 : pr[i] = eos::pressure( r[i]*(e[i] - 0.5*vv) );
1985 [ + + ]: 436651 : if (steady) ma[i] = std::sqrt(vv) / eos::soundspeed( r[i], pr[i] );
1986 : : }
1987 : :
1988 : : std::vector< std::vector< tk::real > > nodefields{
1989 : : std::move(r), std::move(u), std::move(v), std::move(w), std::move(e),
1990 [ - + ][ + + ]: 4179 : std::move(pr) };
[ + - ][ - - ]
[ - - ][ - - ]
1991 [ + + ]: 597 : if (steady) nodefields.push_back( std::move(ma) );
1992 : :
1993 [ + + ]: 695 : for (std::size_t c=0; c<ncomp-5; ++c) {
1994 [ + - ]: 98 : nodefieldnames.push_back( "c" + std::to_string(c) );
1995 [ + - ]: 196 : nodefields.push_back( m_u.extract(5+c) );
1996 : : }
1997 : :
1998 : : // query function to evaluate analytic solution (if defined)
1999 [ + - ]: 597 : auto sol = problems::SOL();
2000 : :
2001 [ + + ]: 597 : if (sol) {
2002 : : const auto& coord = d->Coord();
2003 : : const auto& x = coord[0];
2004 : : const auto& y = coord[1];
2005 : : const auto& z = coord[2];
2006 : : auto an = m_u;
2007 [ + - ][ - - ]: 153 : std::vector< tk::real > ap( m_u.nunk() );
2008 [ + + ]: 15382 : for (std::size_t i=0; i<an.nunk(); ++i) {
2009 [ - + ]: 30458 : auto s = sol( x[i], y[i], z[i], d->T() );
2010 : 15229 : s[1] /= s[0];
2011 : 15229 : s[2] /= s[0];
2012 : 15229 : s[3] /= s[0];
2013 : 15229 : s[4] /= s[0];
2014 [ + + ]: 100476 : for (std::size_t c=0; c<s.size(); ++c) an(i,c) = s[c];
2015 : 15229 : s[4] -= 0.5*(s[1]*s[1] + s[2]*s[2] + s[3]*s[3]);
2016 : 15229 : ap[i] = eos::pressure( s[0]*s[4] );
2017 : : }
2018 [ + + ]: 918 : for (std::size_t c=0; c<5; ++c) {
2019 [ + - ]: 765 : nodefieldnames.push_back( nodefieldnames[c] + "_analytic" );
2020 [ + - ]: 1530 : nodefields.push_back( an.extract(c) );
2021 : : }
2022 [ + - ][ + - ]: 306 : nodefieldnames.push_back( nodefieldnames[5] + "_analytic" );
2023 : : nodefields.push_back( std::move(ap) );
2024 [ + + ]: 251 : for (std::size_t c=0; c<ncomp-5; ++c) {
2025 [ + - ]: 98 : nodefieldnames.push_back( nodefieldnames[6+c] + "_analytic" );
2026 [ + - ][ - - ]: 196 : nodefields.push_back( an.extract(5+c) );
2027 : : }
2028 : : }
2029 : :
2030 : : std::vector< tk::real > bnded, hull, dea;
2031 [ + + ]: 597 : if (g_cfg.get< tag::deactivate >()) {
2032 [ + - ]: 209 : nodefieldnames.push_back( "chbnded" );
2033 [ + - ]: 209 : bnded.resize( m_u.nunk(), 0.0 );
2034 [ + + ]: 1705 : for (const auto& [ch,edges] : m_chbndedge) {
2035 [ + + ]: 794233 : for (const auto& [ed,sint] : edges) {
2036 : 792737 : bnded[ed[0]] = bnded[ed[1]] = 1.0;
2037 : : }
2038 : : }
2039 [ + - ]: 209 : nodefields.push_back( bnded );
2040 : :
2041 [ + - ]: 209 : nodefieldnames.push_back( "activehull" );
2042 [ + - ]: 209 : hull.resize( m_u.nunk(), 0.0 );
2043 [ + + ]: 1705 : for (const auto& [ch,edges] : m_chbndedge) {
2044 : : if (m_inactive.count(ch)) {
2045 [ + + ]: 590585 : for (const auto& [ed,sint] : edges) {
2046 : 589491 : hull[ed[0]] = hull[ed[1]] = 1.0;
2047 : : }
2048 : : }
2049 : : }
2050 [ + - ]: 209 : nodefields.push_back( hull );
2051 : :
2052 [ + - ]: 209 : nodefieldnames.push_back( "deactivated" );
2053 [ + - ][ - - ]: 209 : dea.resize( m_u.nunk(), static_cast<tk::real>(m_deactivated) );
2054 [ + - ]: 209 : nodefields.push_back( dea );
2055 : : }
2056 : :
2057 : : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
2058 : :
2059 : : // Surface output
2060 : :
2061 : : std::vector< std::string > nodesurfnames;
2062 : : std::vector< std::vector< tk::real > > nodesurfs;
2063 : :
2064 : : const auto& f = g_cfg.get< tag::fieldout >();
2065 : :
2066 [ + + ]: 597 : if (!f.empty()) {
2067 [ + - ]: 19 : nodesurfnames.push_back( "density" );
2068 [ + - ]: 19 : nodesurfnames.push_back( "velocityx" );
2069 [ + - ]: 19 : nodesurfnames.push_back( "velocityy" );
2070 [ + - ]: 19 : nodesurfnames.push_back( "velocityz" );
2071 [ + - ]: 19 : nodesurfnames.push_back( "energy" );
2072 [ + - ]: 19 : nodesurfnames.push_back( "pressure" );
2073 : :
2074 [ - + ]: 19 : for (std::size_t c=0; c<ncomp-5; ++c) {
2075 [ - - ]: 0 : nodesurfnames.push_back( "c" + std::to_string(c) );
2076 : : }
2077 : :
2078 [ - + ]: 19 : if (g_cfg.get< tag::deactivate >()) {
2079 [ - - ]: 0 : nodesurfnames.push_back( "chbnded" );
2080 [ - - ]: 0 : nodesurfnames.push_back( "activehull" );
2081 [ - - ]: 0 : nodesurfnames.push_back( "deactivated" );
2082 : : }
2083 : :
2084 [ - + ][ - - ]: 19 : if (steady) nodesurfnames.push_back( "mach" );
2085 : :
2086 [ + - ]: 19 : auto bnode = tk::bfacenodes( m_bface, m_triinpoel );
2087 [ + - ]: 19 : std::set< int > outsets( begin(f), end(f) );
2088 [ + + ]: 71 : for (auto sideset : outsets) {
2089 : : auto b = bnode.find(sideset);
2090 [ + + ]: 52 : if (b == end(bnode)) continue;
2091 : : const auto& nodes = b->second;
2092 : : auto i = nodesurfs.size();
2093 : 46 : auto ns = ncomp + 1;
2094 [ - + ]: 46 : if (g_cfg.get< tag::deactivate >()) ns += 3;
2095 [ - + ]: 46 : if (steady) ++ns;
2096 : : nodesurfs.insert( end(nodesurfs), ns,
2097 [ + - ]: 92 : std::vector< tk::real >( nodes.size() ) );
2098 : : std::size_t j = 0;
2099 [ + + ]: 4876 : for (auto n : nodes) {
2100 [ + - ]: 4830 : const auto s = m_u[n];
2101 : : std::size_t p = 0;
2102 : 4830 : nodesurfs[i+(p++)][j] = s[0];
2103 : 4830 : nodesurfs[i+(p++)][j] = s[1]/s[0];
2104 : 4830 : nodesurfs[i+(p++)][j] = s[2]/s[0];
2105 : 4830 : nodesurfs[i+(p++)][j] = s[3]/s[0];
2106 : 4830 : nodesurfs[i+(p++)][j] = s[4]/s[0];
2107 : 4830 : auto vv = (s[1]*s[1] + s[2]*s[2] + s[3]*s[3])/s[0]/s[0];
2108 : 4830 : auto ei = s[4]/s[0] - 0.5*vv;
2109 : 4830 : auto sp = eos::pressure( s[0]*ei );
2110 : 4830 : nodesurfs[i+(p++)][j] = sp;
2111 [ - + ]: 4830 : for (std::size_t c=0; c<ncomp-5; ++c) nodesurfs[i+(p++)+c][j] = s[5+c];
2112 [ - + ]: 4830 : if (g_cfg.get< tag::deactivate >()) {
2113 : 0 : nodesurfs[i+(p++)][j] = bnded[n];
2114 : 0 : nodesurfs[i+(p++)][j] = hull[n];
2115 : 0 : nodesurfs[i+(p++)][j] = dea[n];
2116 : : }
2117 [ - + ]: 4830 : if (steady) {
2118 : 0 : nodesurfs[i+(p++)][j] = std::sqrt(vv) / eos::soundspeed( s[0], sp );
2119 : : }
2120 : 4830 : ++j;
2121 : : }
2122 : : }
2123 : : }
2124 : :
2125 : : // Send mesh and fields data (solution dump) for output to file
2126 [ + - ][ + - ]: 1194 : d->write( d->Inpoel(), d->Coord(), m_bface, tk::remap(m_bnode,d->Lid()),
2127 [ + - ]: 597 : m_triinpoel, {}, nodefieldnames, {}, nodesurfnames,
2128 : : {}, nodefields, {}, nodesurfs, cb );
2129 [ + - ][ + - ]: 2985 : }
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ][ - - ]
[ - - ]
2130 : :
2131 : : void
2132 : 6040 : ZalCG::out()
2133 : : // *****************************************************************************
2134 : : // Output mesh field data
2135 : : // *****************************************************************************
2136 : : {
2137 : 6040 : auto d = Disc();
2138 : :
2139 : : // Time history
2140 [ + + ][ + + ]: 6040 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ]
2141 : : auto ncomp = m_u.nprop();
2142 : : const auto& inpoel = d->Inpoel();
2143 : 214 : std::vector< std::vector< tk::real > > hist( d->Hist().size() );
2144 : : std::size_t j = 0;
2145 [ + + ]: 282 : for (const auto& p : d->Hist()) {
2146 [ + - ]: 68 : auto e = p.get< tag::elem >(); // host element id
2147 : : const auto& n = p.get< tag::fn >(); // shapefunctions evaluated at point
2148 [ + - ]: 68 : hist[j].resize( ncomp+1, 0.0 );
2149 [ + + ]: 340 : for (std::size_t i=0; i<4; ++i) {
2150 [ + - ]: 272 : const auto u = m_u[ inpoel[e*4+i] ];
2151 : 272 : hist[j][0] += n[i] * u[0];
2152 : 272 : hist[j][1] += n[i] * u[1]/u[0];
2153 : 272 : hist[j][2] += n[i] * u[2]/u[0];
2154 : 272 : hist[j][3] += n[i] * u[3]/u[0];
2155 : 272 : hist[j][4] += n[i] * u[4]/u[0];
2156 : 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];
2157 : 272 : hist[j][5] += n[i] * eos::pressure( u[0]*ei );
2158 [ - + ]: 272 : for (std::size_t c=5; c<ncomp; ++c) hist[j][c+1] += n[i] * u[c];
2159 : : }
2160 : 68 : ++j;
2161 : : }
2162 [ + - ]: 214 : d->history( std::move(hist) );
2163 : 214 : }
2164 : :
2165 : : // Field data
2166 [ + + ][ + + ]: 6040 : if (d->fielditer() or d->fieldtime() or d->fieldrange() or m_finished) {
[ + - ][ + + ]
2167 [ + - ][ + - ]: 2238 : writeFields( CkCallback(CkIndex_ZalCG::integrals(), thisProxy[thisIndex]) );
2168 : : } else {
2169 : 5294 : integrals();
2170 : : }
2171 : 6040 : }
2172 : :
2173 : : void
2174 : 6040 : ZalCG::integrals()
2175 : : // *****************************************************************************
2176 : : // Compute integral quantities for output
2177 : : // *****************************************************************************
2178 : : {
2179 : 6040 : auto d = Disc();
2180 : :
2181 [ + - ][ + - ]: 6040 : if (d->integiter() or d->integtime() or d->integrange()) {
[ - + ]
2182 : :
2183 : : using namespace integrals;
2184 [ - - ]: 0 : std::vector< std::map< int, tk::real > > ints( NUMINT );
2185 : :
2186 : : // Prepend integral vector with metadata on the current time step:
2187 : : // current iteration count, current physical time, time step size
2188 [ - - ][ - - ]: 0 : ints[ ITER ][ 0 ] = static_cast< tk::real >( d->It() );
2189 [ - - ][ - - ]: 0 : ints[ TIME ][ 0 ] = d->T();
2190 [ - - ][ - - ]: 0 : ints[ DT ][ 0 ] = d->Dt();
2191 : :
2192 : : // Compute integrals requested for surfaces requested
2193 : : const auto& reqv = g_cfg.get< tag::integout_integrals >();
2194 : 0 : std::unordered_set< std::string > req( begin(reqv), end(reqv) );
2195 [ - - ][ - - ]: 0 : if (req.count("mass_flow_rate")) {
2196 [ - - ]: 0 : for (const auto& [s,sint] : m_surfint) {
2197 [ - - ]: 0 : auto& mfr = ints[ MASS_FLOW_RATE ][ s ];
2198 : : const auto& nodes = sint.first;
2199 : : const auto& ndA = sint.second;
2200 : : auto n = ndA.data();
2201 [ - - ]: 0 : for (auto p : nodes) {
2202 : 0 : mfr += n[0]*m_u(p,1) + n[1]*m_u(p,2) + n[2]*m_u(p,3);
2203 : 0 : n += 3;
2204 : : }
2205 : : }
2206 : : }
2207 : :
2208 [ - - ]: 0 : auto stream = serialize( d->MeshId(), ints );
2209 [ - - ][ - - ]: 0 : d->contribute( stream.first, stream.second.get(), IntegralsMerger,
2210 [ - - ][ - - ]: 0 : CkCallback(CkIndex_Transporter::integrals(nullptr), d->Tr()) );
2211 : :
2212 : 0 : } else {
2213 : :
2214 : 6040 : step();
2215 : :
2216 : : }
2217 : 6040 : }
2218 : :
2219 : : void
2220 : 5605 : ZalCG::evalLB( int nrestart )
2221 : : // *****************************************************************************
2222 : : // Evaluate whether to do load balancing
2223 : : //! \param[in] nrestart Number of times restarted
2224 : : // *****************************************************************************
2225 : : {
2226 : 5605 : auto d = Disc();
2227 : :
2228 : : // Detect if just returned from a checkpoint and if so, zero timers and
2229 : : // finished flag
2230 [ + + ]: 5605 : if (d->restarted( nrestart )) m_finished = 0;
2231 : :
2232 : : // Load balancing if user frequency is reached or after the second time-step
2233 [ + + ]: 5605 : if (d->lb()) {
2234 : :
2235 : 3614 : AtSync();
2236 [ - + ]: 3614 : if (g_cfg.get< tag::nonblocking >()) dt();
2237 : :
2238 : : } else {
2239 : :
2240 : 1991 : dt();
2241 : :
2242 : : }
2243 : 5605 : }
2244 : :
2245 : : void
2246 : 5600 : ZalCG::evalRestart()
2247 : : // *****************************************************************************
2248 : : // Evaluate whether to save checkpoint/restart
2249 : : // *****************************************************************************
2250 : : {
2251 : 5600 : auto d = Disc();
2252 : :
2253 : 5600 : const auto rsfreq = g_cfg.get< tag::rsfreq >();
2254 : 5600 : const auto benchmark = g_cfg.get< tag::benchmark >();
2255 : :
2256 [ + + ][ - + ]: 5600 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
2257 : :
2258 : 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
2259 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
2260 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
2261 : :
2262 : : } else {
2263 : :
2264 : 5600 : evalLB( /* nrestart = */ -1 );
2265 : :
2266 : : }
2267 : 5600 : }
2268 : :
2269 : : void
2270 : 6040 : ZalCG::step()
2271 : : // *****************************************************************************
2272 : : // Evaluate whether to continue with next time step
2273 : : // *****************************************************************************
2274 : : {
2275 : 6040 : auto d = Disc();
2276 : :
2277 : : // Output one-liner status report to screen
2278 [ + + ]: 6040 : if (thisIndex == 0) d->status();
2279 : :
2280 [ + + ]: 6040 : if (not m_finished) {
2281 : :
2282 : 5600 : evalRestart();
2283 : :
2284 : : } else {
2285 : :
2286 : 440 : auto meshid = d->MeshId();
2287 [ + - ]: 880 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
2288 : 440 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
2289 : :
2290 : : }
2291 : 6040 : }
2292 : :
2293 : : #include "NoWarning/zalcg.def.h"
|
