Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/RieCG.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 RieCG: Riemann, MUSCL, Runge-Kutta, edge-based continuous Galerkin
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "XystBuildConfig.hpp"
14 : : #include "RieCG.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 "Riemann.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 : : //! Runge-Kutta coefficients
41 : : static const std::array< tk::real, 3 > rkcoef{{ 1.0/3.0, 1.0/2.0, 1.0 }};
42 : :
43 : : } // inciter::
44 : :
45 : : using inciter::g_cfg;
46 : : using inciter::RieCG;
47 : :
48 : 559 : RieCG::RieCG( const CProxy_Discretization& disc,
49 : : const std::map< int, std::vector< std::size_t > >& bface,
50 : : const std::map< int, std::vector< std::size_t > >& bnode,
51 : 559 : const std::vector< std::size_t >& triinpoel ) :
52 : : m_disc( disc ),
53 : 559 : m_nrhs( 0 ),
54 : 559 : m_nnorm( 0 ),
55 : 559 : m_nbpint( 0 ),
56 : 559 : m_nbeint( 0 ),
57 : 559 : m_ndeint( 0 ),
58 [ + - ]: 559 : m_ngrad( 0 ),
59 : : m_bnode( bnode ),
60 : : m_bface( bface ),
61 [ + - ][ + - ]: 559 : m_triinpoel( tk::remap( triinpoel, Disc()->Lid() ) ),
62 [ + - ]: 559 : m_u( Disc()->Gid().size(), g_cfg.get< tag::problem_ncomp >() ),
63 : : m_un( m_u.nunk(), m_u.nprop() ),
64 : : m_rhs( m_u.nunk(), m_u.nprop() ),
65 [ + - ]: 559 : m_grad( m_u.nunk(), m_u.nprop()*3 ),
66 : 559 : m_stage( 0 ),
67 [ + - ]: 559 : m_dtp( m_u.nunk(), 0.0 ),
68 [ + - ][ - - ]: 559 : m_tp( m_u.nunk(), g_cfg.get< tag::t0 >() ),
69 [ + - ]: 559 : m_finished( 0 )
70 : : // *****************************************************************************
71 : : // Constructor
72 : : //! \param[in] disc Discretization proxy
73 : : //! \param[in] bface Boundary-faces mapped to side sets used in the input file
74 : : //! \param[in] bnode Boundary-node lists mapped to side sets used in input file
75 : : //! \param[in] triinpoel Boundary-face connectivity where BCs set (global ids)
76 : : // *****************************************************************************
77 : : {
78 : 559 : usesAtSync = true; // enable migration at AtSync
79 : :
80 [ + - ]: 559 : auto d = Disc();
81 : :
82 : : // Create new local ids based on mesh locality
83 : : std::unordered_map< std::size_t, std::size_t > map;
84 : : std::size_t n = 0;
85 : :
86 [ + - ][ + - ]: 559 : auto psup = tk::genPsup( d->Inpoel(), 4, tk::genEsup( d->Inpoel(), 4 ) );
87 [ + + ]: 71217 : for (std::size_t p=0; p<m_u.nunk(); ++p) { // for each point p
88 [ + - ]: 5645 : if (!map.count(p)) map[p] = n++;
89 [ + + ][ + + ]: 802602 : for (auto q : tk::Around(psup,p)) { // for each edge p-q
90 [ + - ]: 65013 : if (!map.count(q)) map[q] = n++;
91 : : }
92 : : }
93 : :
94 : : Assert( map.size() == d->Gid().size(),
95 : : "Mesh-locality reorder map size mismatch" );
96 : :
97 : : // Remap data in bound Discretization object
98 [ + - ]: 559 : d->remap( map );
99 : : // Remap boundary triangle face connectivity
100 [ + - ]: 559 : tk::remap( m_triinpoel, map );
101 : :
102 : : // Compute total box IC volume
103 [ + - ]: 559 : d->boxvol();
104 : :
105 : : // Activate SDAG wait for initially computing integrals
106 [ + - ][ + - ]: 559 : thisProxy[ thisIndex ].wait4int();
107 : 1118 : }
108 : :
109 : : void
110 : 559 : RieCG::setupBC()
111 : : // *****************************************************************************
112 : : // Prepare boundary condition data structures
113 : : // *****************************************************************************
114 : : {
115 : : // Query Dirichlet BC nodes associated to side sets
116 : : std::unordered_map< int, std::unordered_set< std::size_t > > dir;
117 [ + + ]: 1553 : for (const auto& s : g_cfg.get< tag::bc_dir >()) {
118 : : auto k = m_bface.find(s[0]);
119 [ + + ]: 994 : if (k != end(m_bface)) {
120 [ + - ]: 583 : auto& n = dir[ k->first ];
121 [ + + ]: 36287 : for (auto f : k->second) {
122 [ + - ]: 35704 : n.insert( m_triinpoel[f*3+0] );
123 [ + - ]: 35704 : n.insert( m_triinpoel[f*3+1] );
124 [ + - ]: 35704 : n.insert( m_triinpoel[f*3+2] );
125 : : }
126 : : }
127 : : }
128 : :
129 : : // Augment Dirichlet BC nodes with nodes not necessarily part of faces
130 [ + - ]: 559 : const auto& lid = Disc()->Lid();
131 [ + + ]: 1553 : for (const auto& s : g_cfg.get< tag::bc_dir >()) {
132 : : auto k = m_bnode.find(s[0]);
133 [ + + ]: 994 : if (k != end(m_bnode)) {
134 [ + - ]: 597 : auto& n = dir[ k->first ];
135 [ + - ][ + + ]: 37591 : for (auto g : k->second) {
136 : : n.insert( tk::cref_find(lid,g) );
137 : : }
138 : : }
139 : : }
140 : :
141 : : // Collect unique set of nodes + Dirichlet BC components mask
142 : : auto ncomp = m_u.nprop();
143 : 559 : auto nmask = ncomp + 1;
144 : : const auto& dbc = g_cfg.get< tag::bc_dir >();
145 : : std::unordered_map< std::size_t, std::vector< int > > dirbcset;
146 [ + + ]: 1553 : for (const auto& mask : dbc) {
147 [ - + ][ - - ]: 994 : ErrChk( mask.size() == nmask, "Incorrect Dirichlet BC mask ncomp" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
148 : : auto n = dir.find( mask[0] );
149 [ + + ]: 994 : if (n != end(dir))
150 [ + - ][ + + ]: 24739 : for (auto p : n->second) {
151 : : auto& m = dirbcset[p];
152 [ + + ][ + - ]: 24142 : if (m.empty()) m.resize( ncomp, 0 );
153 [ + + ]: 161754 : for (std::size_t c=0; c<ncomp; ++c)
154 [ + + ]: 137612 : if (!m[c]) m[c] = mask[c+1]; // overwrite mask if 0 -> 1
155 : : }
156 : : }
157 : : // Compile streamable list of nodes + Dirichlet BC components mask
158 [ - + ]: 559 : tk::destroy( m_dirbcmasks );
159 [ + + ]: 21055 : for (const auto& [p,mask] : dirbcset) {
160 [ + - ]: 20496 : m_dirbcmasks.push_back( p );
161 [ + - ]: 20496 : m_dirbcmasks.insert( end(m_dirbcmasks), begin(mask), end(mask) );
162 : : }
163 [ - + ][ - - ]: 559 : ErrChk( m_dirbcmasks.size() % nmask == 0, "Dirichlet BC masks incomplete" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
164 : :
165 : : // Query pressure BC nodes associated to side sets
166 : : std::unordered_map< int, std::unordered_set< std::size_t > > pre;
167 [ + + ]: 591 : for (const auto& ss : g_cfg.get< tag::bc_pre >()) {
168 [ + + ]: 64 : for (const auto& s : ss) {
169 : : auto k = m_bface.find(s);
170 [ + + ]: 32 : if (k != end(m_bface)) {
171 [ + - ]: 12 : auto& n = pre[ k->first ];
172 [ + + ]: 132 : for (auto f : k->second) {
173 [ + - ]: 120 : n.insert( m_triinpoel[f*3+0] );
174 [ + - ]: 120 : n.insert( m_triinpoel[f*3+1] );
175 [ + - ]: 120 : n.insert( m_triinpoel[f*3+2] );
176 : : }
177 : : }
178 : : }
179 : : }
180 : :
181 : : // Prepare density and pressure values for pressure BC nodes
182 : : const auto& pbc_set = g_cfg.get< tag::bc_pre >();
183 [ + + ]: 559 : if (!pbc_set.empty()) {
184 : : const auto& pbc_r = g_cfg.get< tag::bc_pre_density >();
185 [ - + ][ - - ]: 16 : ErrChk( pbc_r.size() == pbc_set.size(), "Pressure BC density unspecified" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
186 : : const auto& pbc_p = g_cfg.get< tag::bc_pre_pressure >();
187 [ - + ][ - - ]: 16 : ErrChk( pbc_p.size() == pbc_set.size(), "Pressure BC pressure unspecified" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
188 [ - + ]: 16 : tk::destroy( m_prebcnodes );
189 [ - + ]: 16 : tk::destroy( m_prebcvals );
190 [ + + ]: 28 : for (const auto& [s,n] : pre) {
191 [ + - ]: 12 : m_prebcnodes.insert( end(m_prebcnodes), begin(n), end(n) );
192 [ + + ]: 36 : for (std::size_t p=0; p<pbc_set.size(); ++p) {
193 [ + + ]: 48 : for (auto u : pbc_set[p]) {
194 [ + + ]: 24 : if (s == u) {
195 [ + + ]: 146 : for (std::size_t i=0; i<n.size(); ++i) {
196 [ + - ]: 134 : m_prebcvals.push_back( pbc_r[p] );
197 [ + - ]: 134 : m_prebcvals.push_back( pbc_p[p] );
198 : : }
199 : : }
200 : : }
201 : : }
202 : : }
203 [ - + ][ - - ]: 16 : ErrChk( m_prebcnodes.size()*2 == m_prebcvals.size(),
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
204 : : "Pressure BC data incomplete" );
205 : : }
206 : :
207 : : // Query symmetry BC nodes associated to side sets
208 : : std::unordered_map< int, std::unordered_set< std::size_t > > sym;
209 [ + + ]: 1107 : for (auto s : g_cfg.get< tag::bc_sym >()) {
210 : : auto k = m_bface.find(s);
211 [ + + ]: 548 : if (k != end(m_bface)) {
212 [ + - ]: 337 : auto& n = sym[ k->first ];
213 [ + + ]: 21055 : for (auto f : k->second) {
214 [ + - ]: 20718 : n.insert( m_triinpoel[f*3+0] );
215 [ + - ]: 20718 : n.insert( m_triinpoel[f*3+1] );
216 [ + - ]: 20718 : n.insert( m_triinpoel[f*3+2] );
217 : : }
218 : : }
219 : : }
220 : :
221 : : // Query farfield BC nodes associated to side sets
222 : : std::unordered_map< int, std::unordered_set< std::size_t > > far;
223 [ + + ]: 563 : for (auto s : g_cfg.get< tag::bc_far >()) {
224 : : auto k = m_bface.find(s);
225 [ + + ]: 4 : if (k != end(m_bface)) {
226 [ + - ]: 2 : auto& n = far[ k->first ];
227 [ + + ]: 22 : for (auto f : k->second) {
228 [ + - ]: 20 : n.insert( m_triinpoel[f*3+0] );
229 [ + - ]: 20 : n.insert( m_triinpoel[f*3+1] );
230 [ + - ]: 20 : n.insert( m_triinpoel[f*3+2] );
231 : : }
232 : : }
233 : : }
234 : :
235 : : // Generate unique set of symmetry BC nodes
236 : 559 : tk::destroy( m_symbcnodeset );
237 [ + + ]: 896 : for (const auto& [s,n] : sym) m_symbcnodeset.insert( begin(n), end(n) );
238 : : // Generate unique set of farfield BC nodes
239 : 559 : tk::destroy( m_farbcnodeset );
240 [ + + ]: 561 : for (const auto& [s,n] : far) m_farbcnodeset.insert( begin(n), end(n) );
241 : :
242 : : // If farfield BC is set on a node, will not also set symmetry BC
243 [ + + ]: 583 : for (auto i : m_farbcnodeset) m_symbcnodeset.erase(i);
244 : 559 : }
245 : :
246 : : void
247 : 559 : RieCG::feop()
248 : : // *****************************************************************************
249 : : // Start (re-)computing finite element domain and boundary operators
250 : : // *****************************************************************************
251 : : {
252 : 559 : auto d = Disc();
253 : :
254 : : // Prepare boundary conditions data structures
255 : 559 : setupBC();
256 : :
257 : : // Compute local contributions to boundary normals and integrals
258 : 559 : bndint();
259 : : // Compute local contributions to domain edge integrals
260 : 559 : domint();
261 : :
262 : : // Send boundary point normal contributions to neighbor chares
263 [ + + ]: 559 : if (d->NodeCommMap().empty()) {
264 : 24 : comnorm_complete();
265 : : } else {
266 [ + + ]: 5409 : for (const auto& [c,nodes] : d->NodeCommMap()) {
267 : : decltype(m_bnorm) exp;
268 [ + + ]: 36558 : for (auto i : nodes) {
269 [ + + ]: 103082 : for (const auto& [s,b] : m_bnorm) {
270 : : auto k = b.find(i);
271 [ + + ]: 86451 : if (k != end(b)) exp[s][i] = k->second;
272 : : }
273 : : }
274 [ + - ][ + - ]: 9748 : thisProxy[c].comnorm( exp );
275 : : }
276 : : }
277 : 559 : ownnorm_complete();
278 : 559 : }
279 : :
280 : : void
281 : 559 : RieCG::bndint()
282 : : // *****************************************************************************
283 : : //! Compute local contributions to boundary normals and integrals
284 : : // *****************************************************************************
285 : : {
286 : 559 : auto d = Disc();
287 : : const auto& coord = d->Coord();
288 : : const auto& gid = d->Gid();
289 : : const auto& x = coord[0];
290 : : const auto& y = coord[1];
291 : : const auto& z = coord[2];
292 : :
293 : : // Lambda to compute the inverse distance squared between boundary face
294 : : // centroid and boundary point. Here p is the global node id and c is the
295 : : // the boundary face centroid.
296 : 202176 : auto invdistsq = [&]( const tk::real c[], std::size_t p ){
297 : 202176 : return 1.0 / ( (c[0] - x[p]) * (c[0] - x[p]) +
298 : 202176 : (c[1] - y[p]) * (c[1] - y[p]) +
299 : 202176 : (c[2] - z[p]) * (c[2] - z[p]) );
300 : 559 : };
301 : :
302 : 559 : tk::destroy( m_bnorm );
303 : 559 : tk::destroy( m_bndpoinint );
304 : :
305 [ + + ]: 2021 : for (const auto& [ setid, faceids ] : m_bface) { // for all side sets
306 [ + + ]: 68854 : for (auto f : faceids) { // for all side set triangles
307 : 67392 : const auto N = m_triinpoel.data() + f*3;
308 : : const std::array< tk::real, 3 >
309 : 67392 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
310 : 67392 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
311 : : auto n = tk::cross( ba, ca );
312 : : auto A2 = tk::length( n );
313 : 67392 : n[0] /= A2;
314 : 67392 : n[1] /= A2;
315 : 67392 : n[2] /= A2;
316 : : const tk::real centroid[3] = {
317 : 67392 : (x[N[0]] + x[N[1]] + x[N[2]]) / 3.0,
318 : 67392 : (y[N[0]] + y[N[1]] + y[N[2]]) / 3.0,
319 : 67392 : (z[N[0]] + z[N[1]] + z[N[2]]) / 3.0 };
320 [ + + ]: 269568 : for (const auto& [i,j] : tk::lpoet) {
321 : 202176 : auto p = N[i];
322 : 202176 : tk::real r = invdistsq( centroid, p );
323 : : auto& v = m_bnorm[setid]; // associate side set id
324 : : auto& bpn = v[gid[p]]; // associate global node id of bnd pnt
325 : 202176 : bpn[0] += r * n[0]; // inv.dist.sq-weighted normal
326 : 202176 : bpn[1] += r * n[1];
327 : 202176 : bpn[2] += r * n[2];
328 : 202176 : bpn[3] += r; // inv.dist.sq of node from centroid
329 : : auto& b = m_bndpoinint[gid[p]];// assoc global id of bnd point
330 : 202176 : b[0] += n[0] * A2 / 6.0; // bnd-point integral
331 : 202176 : b[1] += n[1] * A2 / 6.0;
332 : 202176 : b[2] += n[2] * A2 / 6.0;
333 : : }
334 : : }
335 : : }
336 : 559 : }
337 : :
338 : : void
339 : 559 : RieCG::domint()
340 : : // *****************************************************************************
341 : : //! Compute local contributions to domain edge integrals
342 : : // *****************************************************************************
343 : : {
344 : 559 : auto d = Disc();
345 : :
346 : : const auto& gid = d->Gid();
347 : : const auto& inpoel = d->Inpoel();
348 : :
349 : : const auto& coord = d->Coord();
350 : : const auto& x = coord[0];
351 : : const auto& y = coord[1];
352 : : const auto& z = coord[2];
353 : :
354 : 559 : tk::destroy( m_domedgeint );
355 : :
356 [ + + ]: 248768 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
357 : 248209 : const auto N = inpoel.data() + e*4;
358 : : const std::array< tk::real, 3 >
359 : 248209 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
360 : 248209 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
361 : 248209 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
362 : : std::array< std::array< tk::real, 3 >, 4 > grad;
363 : 248209 : grad[1] = tk::cross( ca, da );
364 : 248209 : grad[2] = tk::cross( da, ba );
365 : 248209 : grad[3] = tk::cross( ba, ca );
366 [ + + ]: 992836 : for (std::size_t i=0; i<3; ++i)
367 : 744627 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
368 [ + + ]: 1737463 : for (const auto& [p,q] : tk::lpoed) {
369 [ + + ]: 1489254 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
370 : : tk::real sig = 1.0;
371 [ + + ]: 1489254 : if (ed[0] > ed[1]) {
372 : : std::swap( ed[0], ed[1] );
373 : : sig = -1.0;
374 : : }
375 : : auto& n = m_domedgeint[ ed ];
376 : 1489254 : n[0] += sig * (grad[p][0] - grad[q][0]) / 48.0;
377 : 1489254 : n[1] += sig * (grad[p][1] - grad[q][1]) / 48.0;
378 : 1489254 : n[2] += sig * (grad[p][2] - grad[q][2]) / 48.0;
379 : : }
380 : : }
381 : 559 : }
382 : :
383 : : void
384 : 4874 : RieCG::comnorm( const decltype(m_bnorm)& inbnd )
385 : : // *****************************************************************************
386 : : // Receive contributions to boundary point normals on chare-boundaries
387 : : //! \param[in] inbnd Incoming partial sums of boundary point normals
388 : : // *****************************************************************************
389 : : {
390 : : // Buffer up incoming boundary point normal vector contributions
391 [ + + ]: 9131 : for (const auto& [s,b] : inbnd) {
392 : : auto& bndnorm = m_bnormc[s];
393 [ + + ]: 19310 : for (const auto& [p,n] : b) {
394 : : auto& norm = bndnorm[p];
395 : 15053 : norm[0] += n[0];
396 : 15053 : norm[1] += n[1];
397 : 15053 : norm[2] += n[2];
398 : 15053 : norm[3] += n[3];
399 : : }
400 : : }
401 : :
402 [ + + ]: 4874 : if (++m_nnorm == Disc()->NodeCommMap().size()) {
403 : 535 : m_nnorm = 0;
404 : 535 : comnorm_complete();
405 : : }
406 : 4874 : }
407 : :
408 : : void
409 : 784 : RieCG::registerReducers()
410 : : // *****************************************************************************
411 : : // Configure Charm++ reduction types initiated from this chare array
412 : : //! \details Since this is a [initnode] routine, the runtime system executes the
413 : : //! routine exactly once on every logical node early on in the Charm++ init
414 : : //! sequence. Must be static as it is called without an object. See also:
415 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
416 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
417 : : // *****************************************************************************
418 : : {
419 : 784 : NodeDiagnostics::registerReducers();
420 : 784 : IntegralsMerger = CkReduction::addReducer( integrals::mergeIntegrals );
421 : 784 : }
422 : :
423 : : void
424 : : // cppcheck-suppress unusedFunction
425 : 9039 : RieCG::ResumeFromSync()
426 : : // *****************************************************************************
427 : : // Return from migration
428 : : //! \details This is called when load balancing (LB) completes. The presence of
429 : : //! this function does not affect whether or not we block on LB.
430 : : // *****************************************************************************
431 : : {
432 [ - + ][ - - ]: 9039 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
433 : :
434 [ + - ]: 9039 : if (!g_cfg.get< tag::nonblocking >()) dt();
435 : 9039 : }
436 : :
437 : : void
438 : 559 : RieCG::setup( tk::real v )
439 : : // *****************************************************************************
440 : : // Start setup for solution
441 : : //! \param[in] v Total volume within user-specified box
442 : : // *****************************************************************************
443 : : {
444 : 559 : auto d = Disc();
445 : :
446 : : // Store user-defined box IC volume
447 : 559 : Disc()->Boxvol() = v;
448 : :
449 : : // Set initial conditions
450 : 559 : problems::initialize( d->Coord(), m_u, d->T(), d->BoxNodes() );
451 : :
452 : : // Query time history field output labels from all PDEs integrated
453 [ + + ]: 559 : if (!g_cfg.get< tag::histout >().empty()) {
454 : : std::vector< std::string > var
455 [ - + ][ + + ]: 210 : {"density", "xvelocity", "yvelocity", "zvelocity", "energy", "pressure"};
[ - + ][ - - ]
[ - - ]
456 : : auto ncomp = m_u.nprop();
457 [ + + ]: 34 : for (std::size_t c=5; c<ncomp; ++c)
458 [ + - ]: 8 : var.push_back( "c" + std::to_string(c-5) );
459 [ + - ]: 30 : d->histheader( std::move(var) );
460 : 30 : }
461 : :
462 : : // Compute finite element operators
463 : 559 : feop();
464 [ + - ][ + - ]: 619 : }
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
465 : :
466 : : void
467 : 559 : RieCG::start()
468 : : // *****************************************************************************
469 : : // Start time stepping
470 : : // *****************************************************************************
471 : : {
472 : : // Set flag that indicates that we are now during time stepping
473 : 559 : Disc()->Initial( 0 );
474 : : // Start timer measuring time stepping wall clock time
475 : 559 : Disc()->Timer().zero();
476 : : // Zero grind-timer
477 : 559 : Disc()->grindZero();
478 : : // Continue to first time step
479 : 559 : dt();
480 : 559 : }
481 : :
482 : : void
483 : 559 : RieCG::bnorm()
484 : : // *****************************************************************************
485 : : // Combine own and communicated portions of the boundary point normals
486 : : // *****************************************************************************
487 : : {
488 : 559 : const auto& lid = Disc()->Lid();
489 : :
490 : : // Combine own and communicated contributions to boundary point normals
491 [ + + ]: 1913 : for (const auto& [s,b] : m_bnormc) {
492 : : auto& bndnorm = m_bnorm[s];
493 [ + + ]: 13104 : for (const auto& [g,n] : b) {
494 : : auto& norm = bndnorm[g];
495 : 11750 : norm[0] += n[0];
496 : 11750 : norm[1] += n[1];
497 : 11750 : norm[2] += n[2];
498 : 11750 : norm[3] += n[3];
499 : : }
500 : : }
501 : 559 : tk::destroy( m_bnormc );
502 : :
503 : : // Divide summed point normals by the sum of the inverse distance squared
504 [ + + ]: 2089 : for (auto& [s,b] : m_bnorm) {
505 [ + + ]: 49173 : for (auto& [g,n] : b) {
506 : 47643 : n[0] /= n[3];
507 : 47643 : n[1] /= n[3];
508 : 47643 : n[2] /= n[3];
509 : : Assert( (n[0]*n[0] + n[1]*n[1] + n[2]*n[2] - 1.0) <
510 : : 1.0e+3*std::numeric_limits< tk::real >::epsilon(),
511 : : "Non-unit normal" );
512 : : }
513 : : }
514 : :
515 : : // Replace global->local ids associated to boundary point normals
516 : : decltype(m_bnorm) loc;
517 [ + + ]: 2089 : for (auto& [s,b] : m_bnorm) {
518 : : auto& bnd = loc[s];
519 [ + + ]: 49173 : for (auto&& [g,n] : b) {
520 : 47643 : bnd[ tk::cref_find(lid,g) ] = std::move(n);
521 : : }
522 : : }
523 : : m_bnorm = std::move(loc);
524 : 559 : }
525 : :
526 : : void
527 : 559 : RieCG::streamable()
528 : : // *****************************************************************************
529 : : // Convert integrals into streamable data structures
530 : : // *****************************************************************************
531 : : {
532 : : // Generate boundary element symmetry BC flags
533 : 559 : m_besym.resize( m_triinpoel.size() );
534 : : std::size_t i = 0;
535 [ + + ]: 202735 : for (auto p : m_triinpoel) {
536 : 202176 : m_besym[i++] = static_cast< std::uint8_t >(m_symbcnodeset.count(p));
537 : : }
538 : :
539 : : // Query surface integral output nodes
540 : : std::unordered_map< int, std::vector< std::size_t > > surfintnodes;
541 : : const auto& is = g_cfg.get< tag::integout >();
542 [ + - ]: 559 : std::set< int > outsets( begin(is), end(is) );
543 [ + + ]: 587 : for (auto s : outsets) {
544 : : auto m = m_bface.find(s);
545 [ + + ]: 28 : if (m != end(m_bface)) {
546 [ + - ]: 10 : auto& n = surfintnodes[ m->first ]; // associate set id
547 [ + + ]: 110 : for (auto f : m->second) { // face ids on side set
548 [ + - ]: 100 : n.push_back( m_triinpoel[f*3+0] ); // nodes on side set
549 [ + - ]: 100 : n.push_back( m_triinpoel[f*3+1] );
550 [ + - ]: 100 : n.push_back( m_triinpoel[f*3+2] );
551 : : }
552 : : }
553 : : }
554 [ + - ][ + + ]: 569 : for (auto& [s,n] : surfintnodes) tk::unique( n );
555 : : // Prepare surface integral data
556 : 559 : tk::destroy( m_surfint );
557 [ + - ]: 559 : const auto& gid = Disc()->Gid();
558 [ + + ]: 569 : for (auto&& [s,n] : surfintnodes) {
559 [ + - ]: 10 : auto& sint = m_surfint[s]; // associate set id
560 : 10 : auto& nodes = sint.first;
561 : 10 : auto& ndA = sint.second;
562 : : nodes = std::move(n);
563 [ + - ]: 10 : ndA.resize( nodes.size()*3 );
564 : : std::size_t a = 0;
565 [ + + ]: 120 : for (auto p : nodes) {
566 : : const auto& b = tk::cref_find( m_bndpoinint, gid[p] );
567 : 110 : ndA[a*3+0] = b[0]; // store ni * dA
568 : 110 : ndA[a*3+1] = b[1];
569 : 110 : ndA[a*3+2] = b[2];
570 : 110 : ++a;
571 : : }
572 : : }
573 : 559 : tk::destroy( m_bndpoinint );
574 : :
575 : : // Generate domain superedges
576 [ + - ]: 559 : domsuped();
577 : 559 : tk::destroy( m_domedgeint );
578 : :
579 : : // Convert symmetry BC data to streamable data structures
580 [ - + ]: 559 : tk::destroy( m_symbcnodes );
581 [ - + ]: 559 : tk::destroy( m_symbcnorms );
582 [ + + ]: 12606 : for (auto p : m_symbcnodeset) {
583 [ + + ]: 54047 : for (const auto& s : g_cfg.get< tag::bc_sym >()) {
584 : : auto m = m_bnorm.find(s);
585 [ + + ]: 42000 : if (m != end(m_bnorm)) {
586 : : auto r = m->second.find(p);
587 [ + + ]: 39485 : if (r != end(m->second)) {
588 [ + - ]: 14446 : m_symbcnodes.push_back( p );
589 [ + - ]: 14446 : m_symbcnorms.push_back( r->second[0] );
590 [ + - ]: 14446 : m_symbcnorms.push_back( r->second[1] );
591 [ + - ]: 14446 : m_symbcnorms.push_back( r->second[2] );
592 : : }
593 : : }
594 : : }
595 : : }
596 : 559 : tk::destroy( m_symbcnodeset );
597 : :
598 : : // Convert farfield BC data to streamable data structures
599 [ - + ]: 559 : tk::destroy( m_farbcnodes );
600 [ - + ]: 559 : tk::destroy( m_farbcnorms );
601 [ + + ]: 583 : for (auto p : m_farbcnodeset) {
602 [ + + ]: 72 : for (const auto& s : g_cfg.get< tag::bc_far >()) {
603 : : auto n = m_bnorm.find(s);
604 [ + + ]: 48 : if (n != end(m_bnorm)) {
605 : : auto a = n->second.find(p);
606 [ + - ]: 24 : if (a != end(n->second)) {
607 [ + - ]: 24 : m_farbcnodes.push_back( p );
608 [ + - ]: 24 : m_farbcnorms.push_back( a->second[0] );
609 [ + - ]: 24 : m_farbcnorms.push_back( a->second[1] );
610 [ + - ]: 24 : m_farbcnorms.push_back( a->second[2] );
611 : : }
612 : : }
613 : : }
614 : : }
615 : 559 : tk::destroy( m_farbcnodeset );
616 : 559 : tk::destroy( m_bnorm );
617 : 559 : }
618 : :
619 : : void
620 : 559 : RieCG::domsuped()
621 : : // *****************************************************************************
622 : : // Generate superedge-groups for domain-edge loops
623 : : //! \see See Lohner, Sec. 15.1.6.2, An Introduction to Applied CFD Techniques,
624 : : //! Wiley, 2008.
625 : : // *****************************************************************************
626 : : {
627 : : Assert( !m_domedgeint.empty(), "No domain edges to group" );
628 : :
629 : : #ifndef NDEBUG
630 : : auto nedge = m_domedgeint.size();
631 : : #endif
632 : :
633 : 559 : const auto& inpoel = Disc()->Inpoel();
634 : 559 : const auto& lid = Disc()->Lid();
635 : 559 : const auto& gid = Disc()->Gid();
636 : :
637 : : tk::destroy( m_dsupedge[0] );
638 : : tk::destroy( m_dsupedge[1] );
639 : : tk::destroy( m_dsupedge[2] );
640 : :
641 : : tk::destroy( m_dsupint[0] );
642 : : tk::destroy( m_dsupint[1] );
643 : : tk::destroy( m_dsupint[2] );
644 : :
645 : : tk::UnsMesh::FaceSet untri;
646 [ + + ]: 248768 : for (std::size_t e=0; e<inpoel.size()/4; e++) {
647 : : std::size_t N[4] = {
648 : 248209 : inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
649 [ + - ][ + + ]: 1241045 : for (const auto& [a,b,c] : tk::lpofa) untri.insert( { N[a], N[b], N[c] } );
650 : : }
651 : :
652 [ + + ]: 248768 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
653 : : std::size_t N[4] = {
654 : 248209 : inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
655 : : int f = 0;
656 : : tk::real sig[6];
657 [ + + ]: 1737463 : decltype(m_domedgeint)::const_iterator d[6];
658 [ + + ]: 710288 : for (const auto& [p,q] : tk::lpoed) {
659 [ + + ]: 676370 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
660 [ + + ]: 961133 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
661 : 676370 : d[f] = m_domedgeint.find( ed );
662 [ + + ]: 676370 : if (d[f] == end(m_domedgeint)) break; else ++f;
663 : : }
664 [ + + ]: 248209 : if (f == 6) {
665 [ + - ]: 33918 : m_dsupedge[0].push_back( N[0] );
666 [ + - ]: 33918 : m_dsupedge[0].push_back( N[1] );
667 [ + - ]: 33918 : m_dsupedge[0].push_back( N[2] );
668 [ + - ]: 33918 : m_dsupedge[0].push_back( N[3] );
669 [ + + ]: 169590 : for (const auto& [a,b,c] : tk::lpofa) untri.erase( { N[a], N[b], N[c] } );
670 [ + + ]: 237426 : for (int ed=0; ed<6; ++ed) {
671 [ + - ]: 203508 : m_dsupint[0].push_back( sig[ed] * d[ed]->second[0] );
672 [ + - ]: 203508 : m_dsupint[0].push_back( sig[ed] * d[ed]->second[1] );
673 [ + - ]: 203508 : m_dsupint[0].push_back( sig[ed] * d[ed]->second[2] );
674 : 203508 : m_domedgeint.erase( d[ed] );
675 : : }
676 : : }
677 : : }
678 : :
679 [ + + ]: 408969 : for (const auto& N : untri) {
680 : : int f = 0;
681 : : tk::real sig[3];
682 [ + + ]: 1633640 : decltype(m_domedgeint)::const_iterator d[3];
683 [ + + ]: 669667 : for (const auto& [p,q] : tk::lpoet) {
684 [ + + ]: 645191 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
685 [ + + ]: 986590 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
686 : 645191 : d[f] = m_domedgeint.find( ed );
687 [ + + ]: 645191 : if (d[f] == end(m_domedgeint)) break; else ++f;
688 : : }
689 [ + + ]: 408410 : if (f == 3) {
690 [ + - ]: 24476 : m_dsupedge[1].push_back( N[0] );
691 [ + - ]: 24476 : m_dsupedge[1].push_back( N[1] );
692 [ + - ]: 24476 : m_dsupedge[1].push_back( N[2] );
693 [ + + ]: 97904 : for (int ed=0; ed<3; ++ed) {
694 [ + - ]: 73428 : m_dsupint[1].push_back( sig[ed] * d[ed]->second[0] );
695 [ + - ]: 73428 : m_dsupint[1].push_back( sig[ed] * d[ed]->second[1] );
696 [ + - ]: 73428 : m_dsupint[1].push_back( sig[ed] * d[ed]->second[2] );
697 : 73428 : m_domedgeint.erase( d[ed] );
698 : : }
699 : : }
700 : : }
701 : :
702 [ + - ]: 559 : m_dsupedge[2].resize( m_domedgeint.size()*2 );
703 [ + - ]: 559 : m_dsupint[2].resize( m_domedgeint.size()*3 );
704 : : std::size_t k = 0;
705 [ + + ]: 89595 : for (const auto& [ed,d] : m_domedgeint) {
706 : 89036 : auto e = m_dsupedge[2].data() + k*2;
707 : 89036 : e[0] = tk::cref_find( lid, ed[0] );
708 : 89036 : e[1] = tk::cref_find( lid, ed[1] );
709 : 89036 : auto i = m_dsupint[2].data() + k*3;
710 : 89036 : i[0] = d[0];
711 : 89036 : i[1] = d[1];
712 : 89036 : i[2] = d[2];
713 : 89036 : ++k;
714 : : }
715 : :
716 : : //std::cout << std::setprecision(2)
717 : : // << "superedges: ntet:" << m_dsupedge[0].size()/4 << "(nedge:"
718 : : // << m_dsupedge[0].size()/4*6 << ","
719 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[0].size()/4*6 ) /
720 : : // static_cast< tk::real >( nedge )
721 : : // << "%) + ntri:" << m_dsupedge[1].size()/3
722 : : // << "(nedge:" << m_dsupedge[1].size() << ","
723 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[1].size() ) /
724 : : // static_cast< tk::real >( nedge )
725 : : // << "%) + nedge:"
726 : : // << m_dsupedge[2].size()/2 << "("
727 : : // << 100.0 * static_cast< tk::real >( m_dsupedge[2].size()/2 ) /
728 : : // static_cast< tk::real >( nedge )
729 : : // << "%) = " << m_dsupedge[0].size()/4*6 + m_dsupedge[1].size() +
730 : : // m_dsupedge[2].size()/2 << " of "<< nedge << " total edges\n";
731 : :
732 : : Assert( m_dsupedge[0].size()/4*6 + m_dsupedge[1].size() +
733 : : m_dsupedge[2].size()/2 == nedge,
734 : : "Not all edges accounted for in superedge groups" );
735 : 559 : }
736 : :
737 : : void
738 : : // cppcheck-suppress unusedFunction
739 : 559 : RieCG::merge()
740 : : // *****************************************************************************
741 : : // Combine own and communicated portions of the integrals
742 : : // *****************************************************************************
743 : : {
744 : 559 : auto d = Disc();
745 : :
746 : : // Combine own and communicated contributions to boundary point normals
747 : 559 : bnorm();
748 : :
749 : : // Convert integrals into streamable data structures
750 : 559 : streamable();
751 : :
752 : : // Enforce boundary conditions using (re-)computed boundary data
753 : 559 : BC( d->T() );
754 : :
755 [ + - ]: 559 : if (d->Initial()) {
756 : : // Output initial conditions to file
757 [ + - ][ + - ]: 1677 : writeFields( CkCallback(CkIndex_RieCG::start(), thisProxy[thisIndex]) );
758 : : } else {
759 : 0 : feop_complete();
760 : : }
761 : 559 : }
762 : :
763 : : void
764 : 38599 : RieCG::BC( tk::real t )
765 : : // *****************************************************************************
766 : : // Apply boundary conditions
767 : : //! \param[in] t Physical time
768 : : // *****************************************************************************
769 : : {
770 : 38599 : auto d = Disc();
771 : :
772 : : // Apply Dirichlet BCs
773 [ + - ]: 38599 : physics::dirbc( m_u, t, d->Coord(), d->BoxNodes(), m_dirbcmasks );
774 : :
775 : : // Apply symmetry BCs
776 : 38599 : physics::symbc( m_u, m_symbcnodes, m_symbcnorms, /*pos=*/1 );
777 : :
778 : : // Apply farfield BCs
779 : 38599 : physics::farbc( m_u, m_farbcnodes, m_farbcnorms );
780 : :
781 : : // Apply pressure BCs
782 : 38599 : physics::prebc( m_u, m_prebcnodes, m_prebcvals );
783 : 38599 : }
784 : :
785 : : void
786 : 12680 : RieCG::dt()
787 : : // *****************************************************************************
788 : : // Compute time step size
789 : : // *****************************************************************************
790 : : {
791 : 12680 : tk::real mindt = std::numeric_limits< tk::real >::max();
792 : :
793 : 12680 : auto const_dt = g_cfg.get< tag::dt >();
794 : : auto eps = std::numeric_limits< tk::real >::epsilon();
795 : 12680 : auto d = Disc();
796 : :
797 : : // use constant dt if configured
798 [ + + ]: 12680 : if (std::abs(const_dt) > eps) {
799 : :
800 : : // cppcheck-suppress redundantInitialization
801 : 2920 : mindt = const_dt;
802 : :
803 : : } else {
804 : :
805 : : const auto& vol = d->Vol();
806 : 9760 : auto cfl = g_cfg.get< tag::cfl >();
807 : :
808 [ + + ]: 9760 : if (g_cfg.get< tag::steady >()) {
809 : :
810 [ + + ]: 5620 : for (std::size_t p=0; p<m_u.nunk(); ++p) {
811 : 5570 : auto r = m_u(p,0);
812 : 5570 : auto u = m_u(p,1)/r;
813 : 5570 : auto v = m_u(p,2)/r;
814 : 5570 : auto w = m_u(p,3)/r;
815 [ - + ]: 5570 : auto pr = eos::pressure( m_u(p,4) - 0.5*r*(u*u + v*v + w*w) );
816 [ - + ][ - + ]: 5570 : auto c = eos::soundspeed( r, std::max(pr,0.0) );
817 : 5570 : auto L = std::cbrt( vol[p] );
818 : 5570 : auto vel = std::sqrt( u*u + v*v + w*w );
819 [ - + ]: 5570 : m_dtp[p] = L / std::max( vel+c, 1.0e-8 ) * cfl;
820 : : }
821 : 50 : mindt = *std::min_element( begin(m_dtp), end(m_dtp) );
822 : :
823 : : } else {
824 : :
825 [ + + ]: 1119315 : for (std::size_t p=0; p<m_u.nunk(); ++p) {
826 : 1109605 : auto r = m_u(p,0);
827 : 1109605 : auto u = m_u(p,1)/r;
828 : 1109605 : auto v = m_u(p,2)/r;
829 : 1109605 : auto w = m_u(p,3)/r;
830 [ - + ]: 1109605 : auto pr = eos::pressure( m_u(p,4) - 0.5*r*(u*u + v*v + w*w) );
831 [ - + ][ - + ]: 1109605 : auto c = eos::soundspeed( r, std::max(pr,0.0) );
832 : 1109605 : auto L = std::cbrt( vol[p] );
833 : 1109605 : auto vel = std::sqrt( u*u + v*v + w*w );
834 [ - + ][ + + ]: 1109605 : auto euler_dt = L / std::max( vel+c, 1.0e-8 );
835 : 1109605 : mindt = std::min( mindt, euler_dt );
836 : : }
837 : 9710 : mindt *= cfl;
838 : :
839 : : }
840 : :
841 : : }
842 : :
843 : : // Actiavate SDAG waits for next time step stage
844 [ + - ]: 12680 : thisProxy[ thisIndex ].wait4grad();
845 [ + - ]: 12680 : thisProxy[ thisIndex ].wait4rhs();
846 : :
847 : : // Contribute to minimum dt across all chares and advance to next step
848 [ + - ]: 12680 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
849 : 12680 : CkCallback(CkReductionTarget(RieCG,advance), thisProxy) );
850 : 12680 : }
851 : :
852 : : void
853 : 12680 : RieCG::advance( tk::real newdt )
854 : : // *****************************************************************************
855 : : // Advance equations to next time step
856 : : //! \param[in] newdt The smallest dt across the whole problem
857 : : // *****************************************************************************
858 : : {
859 : : // Detect blowup
860 : : auto eps = std::numeric_limits< tk::real >::epsilon();
861 [ - + ]: 12680 : if (newdt < eps) m_finished = 1;
862 : :
863 : : // Set new time step size
864 [ + - ]: 12680 : if (m_stage == 0) Disc()->setdt( newdt );
865 : :
866 : 12680 : grad();
867 : 12680 : }
868 : :
869 : : void
870 : 38040 : RieCG::grad()
871 : : // *****************************************************************************
872 : : // Compute gradients for next time step
873 : : // *****************************************************************************
874 : : {
875 : 38040 : auto d = Disc();
876 : :
877 : 38040 : riemann::grad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_grad );
878 : :
879 : : // Send gradient contributions to neighbor chares
880 [ + + ]: 38040 : if (d->NodeCommMap().empty()) {
881 : 1404 : comgrad_complete();
882 : : } else {
883 : : const auto& lid = d->Lid();
884 [ + + ]: 332910 : for (const auto& [c,n] : d->NodeCommMap()) {
885 : : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
886 [ + - ][ + + ]: 2119578 : for (auto g : n) exp[g] = m_grad[ tk::cref_find(lid,g) ];
887 [ + - ][ + - ]: 592548 : thisProxy[c].comgrad( exp );
888 : : }
889 : : }
890 : 38040 : owngrad_complete();
891 : 38040 : }
892 : :
893 : : void
894 : 296274 : RieCG::comgrad(
895 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& ingrad )
896 : : // *****************************************************************************
897 : : // Receive contributions to node gradients on chare-boundaries
898 : : //! \param[in] ingrad Partial contributions to chare-boundary nodes. Key:
899 : : //! global mesh node IDs, value: contributions for all scalar components.
900 : : //! \details This function receives contributions to m_grad, which stores the
901 : : //! gradients at mesh nodes. While m_grad stores own contributions, m_gradc
902 : : //! collects the neighbor chare contributions during communication. This way
903 : : //! work on m_grad and m_gradc is overlapped. The two are combined in rhs().
904 : : // *****************************************************************************
905 : : {
906 : : using tk::operator+=;
907 [ + + ]: 3942882 : for (const auto& [g,r] : ingrad) m_gradc[g] += r;
908 : :
909 : : // When we have heard from all chares we communicate with, this chare is done
910 [ + + ]: 296274 : if (++m_ngrad == Disc()->NodeCommMap().size()) {
911 : 36636 : m_ngrad = 0;
912 : 36636 : comgrad_complete();
913 : : }
914 : 296274 : }
915 : :
916 : : void
917 : 38040 : RieCG::rhs()
918 : : // *****************************************************************************
919 : : // Compute right-hand side of transport equations
920 : : // *****************************************************************************
921 : : {
922 : 38040 : auto d = Disc();
923 : : const auto& lid = d->Lid();
924 : 38040 : const auto steady = g_cfg.get< tag::steady >();
925 : :
926 : : // Combine own and communicated contributions to gradients
927 [ + + ]: 1127304 : for (const auto& [g,r] : m_gradc) {
928 : 1089264 : auto i = tk::cref_find( lid, g );
929 [ + + ]: 17931342 : for (std::size_t c=0; c<r.size(); ++c) m_grad(i,c) += r[c];
930 : : }
931 : 38040 : tk::destroy(m_gradc);
932 : :
933 : : // divide weak result in gradients by nodal volume
934 : : const auto& vol = d->Vol();
935 [ + + ]: 3614505 : for (std::size_t p=0; p<m_grad.nunk(); ++p)
936 [ + + ]: 60027930 : for (std::size_t c=0; c<m_grad.nprop(); ++c)
937 : 56451465 : m_grad(p,c) /= vol[p];
938 : :
939 : : // Compute own portion of right-hand side for all equations
940 [ + + ]: 38040 : auto prev_rkcoef = m_stage == 0 ? 0.0 : rkcoef[m_stage-1];
941 : :
942 [ + + ]: 38040 : if (steady) {
943 [ + + ]: 16860 : for (std::size_t p=0; p<m_tp.size(); ++p) m_tp[p] += prev_rkcoef * m_dtp[p];
944 : : }
945 : :
946 : 38040 : riemann::rhs( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_besym, m_grad,
947 : 38040 : m_u, d->V(), d->T(), m_tp, m_rhs );
948 : :
949 [ + + ]: 38040 : if (steady) {
950 [ + + ]: 16860 : for (std::size_t p=0; p<m_tp.size(); ++p) m_tp[p] -= prev_rkcoef * m_dtp[p];
951 : : }
952 : :
953 : : // Communicate rhs to other chares on chare-boundary
954 [ + + ]: 38040 : if (d->NodeCommMap().empty()) {
955 : 1404 : comrhs_complete();
956 : : } else {
957 [ + + ]: 332910 : for (const auto& [c,n] : d->NodeCommMap()) {
958 : : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
959 [ + - ][ + + ]: 2119578 : for (auto g : n) exp[g] = m_rhs[ tk::cref_find(lid,g) ];
960 [ + - ][ + - ]: 592548 : thisProxy[c].comrhs( exp );
961 : : }
962 : : }
963 : 38040 : ownrhs_complete();
964 : 38040 : }
965 : :
966 : : void
967 : 296274 : RieCG::comrhs(
968 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& inrhs )
969 : : // *****************************************************************************
970 : : // Receive contributions to right-hand side vector on chare-boundaries
971 : : //! \param[in] inrhs Partial contributions of RHS to chare-boundary nodes. Key:
972 : : //! global mesh node IDs, value: contributions for all scalar components.
973 : : //! \details This function receives contributions to m_rhs, which stores the
974 : : //! right hand side vector at mesh nodes. While m_rhs stores own
975 : : //! contributions, m_rhsc collects the neighbor chare contributions during
976 : : //! communication. This way work on m_rhs and m_rhsc is overlapped. The two
977 : : //! are combined in solve().
978 : : // *****************************************************************************
979 : : {
980 : : using tk::operator+=;
981 [ + + ]: 3942882 : for (const auto& [g,r] : inrhs) m_rhsc[g] += r;
982 : :
983 : : // When we have heard from all chares we communicate with, this chare is done
984 [ + + ]: 296274 : if (++m_nrhs == Disc()->NodeCommMap().size()) {
985 : 36636 : m_nrhs = 0;
986 : 36636 : comrhs_complete();
987 : : }
988 : 296274 : }
989 : :
990 : : void
991 : : // cppcheck-suppress unusedFunction
992 : 38040 : RieCG::solve()
993 : : // *****************************************************************************
994 : : // Advance systems of equations
995 : : // *****************************************************************************
996 : : {
997 : 38040 : auto d = Disc();
998 : : const auto lid = d->Lid();
999 : 38040 : const auto steady = g_cfg.get< tag::steady >();
1000 : :
1001 : : // Combine own and communicated contributions to rhs
1002 [ + + ]: 1127304 : for (const auto& [g,r] : m_rhsc) {
1003 : 1089264 : auto i = tk::cref_find( lid, g );
1004 [ + + ]: 6703290 : for (std::size_t c=0; c<r.size(); ++c) m_rhs(i,c) += r[c];
1005 : : }
1006 : 38040 : tk::destroy(m_rhsc);
1007 : :
1008 : : // Update state at time n
1009 [ + + ]: 38040 : if (m_stage == 0) m_un = m_u;
1010 : :
1011 : : // Advance solution
1012 : : auto dt = d->Dt();
1013 : : const auto& vol = d->Vol();
1014 [ + + ]: 3614505 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
1015 [ + + ]: 3576465 : if (steady) dt = m_dtp[i];
1016 [ + + ]: 22393620 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
1017 : 18817155 : m_u(i,c) = m_un(i,c) - rkcoef[m_stage] * dt * m_rhs(i,c) / vol[i];
1018 : : }
1019 : : }
1020 : :
1021 : : // Configure and apply scalar source to solution (if defined)
1022 [ + - ]: 38040 : auto src = problems::PHYS_SRC();
1023 [ + + ][ + - ]: 38040 : if (src) src( d->Coord(), d->T(), m_u );
1024 : :
1025 : : // Enforce boundary conditions
1026 [ + - ]: 38040 : BC( d->T() + rkcoef[m_stage] * d->Dt() );
1027 : :
1028 [ + + ]: 38040 : if (m_stage < 2) {
1029 : :
1030 : : // Activate SDAG wait for next time step stage
1031 [ + - ][ + - ]: 25360 : thisProxy[ thisIndex ].wait4grad();
1032 [ + - ][ + - ]: 25360 : thisProxy[ thisIndex ].wait4rhs();
1033 : :
1034 : : // start next time step stage
1035 [ + - ]: 25360 : stage();
1036 : :
1037 : : } else {
1038 : :
1039 : : // Activate SDAG waits for finishing this time step stage
1040 [ + - ][ + - ]: 12680 : thisProxy[ thisIndex ].wait4stage();
1041 : : // Compute diagnostics, e.g., residuals
1042 : 12680 : auto diag_iter = g_cfg.get< tag::diag_iter >();
1043 [ + - ]: 12680 : auto diag = m_diag.rhocompute( *d, m_u, m_un, diag_iter );
1044 : : // Increase number of iterations and physical time
1045 [ + - ]: 12680 : d->next();
1046 : : // Advance physical time for local time stepping
1047 [ + + ]: 12680 : if (steady) {
1048 : : using tk::operator+=;
1049 [ + - ]: 50 : m_tp += m_dtp;
1050 : : }
1051 : : // Evaluate residuals
1052 [ + + ][ + - ]: 13302 : if (!diag) evalres( std::vector< tk::real >( m_u.nprop(), 1.0 ) );
[ + - ][ - - ]
1053 : :
1054 : : }
1055 : 38040 : }
1056 : :
1057 : : void
1058 : 12680 : RieCG::evalres( const std::vector< tk::real >& l2res )
1059 : : // *****************************************************************************
1060 : : // Evaluate residuals
1061 : : //! \param[in] l2res L2-norms of the residual for each scalar component
1062 : : //! computed across the whole problem
1063 : : // *****************************************************************************
1064 : : {
1065 [ + + ]: 12680 : if (g_cfg.get< tag::steady >()) {
1066 : 50 : const auto rc = g_cfg.get< tag::rescomp >() - 1;
1067 : 50 : Disc()->residual( l2res[rc] );
1068 : : }
1069 : :
1070 : 12680 : refine();
1071 : 12680 : }
1072 : :
1073 : : void
1074 : 12680 : RieCG::refine()
1075 : : // *****************************************************************************
1076 : : // Optionally refine/derefine mesh
1077 : : // *****************************************************************************
1078 : : {
1079 : 12680 : auto d = Disc();
1080 : :
1081 : : // See if this is the last time step
1082 [ + + ]: 12680 : if (d->finished()) m_finished = 1;
1083 : :
1084 : 12680 : auto dtref = g_cfg.get< tag::href_dt >();
1085 : 12680 : auto dtfreq = g_cfg.get< tag::href_dtfreq >();
1086 : :
1087 : : // if t>0 refinement enabled and we hit the frequency
1088 [ - + ][ - - ]: 12680 : if (dtref && !(d->It() % dtfreq)) { // refine
1089 : :
1090 : 0 : d->refined() = 1;
1091 : 0 : d->startvol();
1092 : 0 : d->Ref()->dtref( m_bface, m_bnode, m_triinpoel );
1093 : :
1094 : : // Activate SDAG waits for re-computing the integrals
1095 [ - - ]: 0 : thisProxy[ thisIndex ].wait4int();
1096 : :
1097 : : } else { // do not refine
1098 : :
1099 : 12680 : d->refined() = 0;
1100 : 12680 : feop_complete();
1101 : 12680 : resize_complete();
1102 : :
1103 : : }
1104 : 12680 : }
1105 : :
1106 : : void
1107 : 0 : RieCG::resizePostAMR(
1108 : : const std::vector< std::size_t >& /*ginpoel*/,
1109 : : const tk::UnsMesh::Chunk& chunk,
1110 : : const tk::UnsMesh::Coords& coord,
1111 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
1112 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
1113 : : const std::set< std::size_t >& removedNodes,
1114 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
1115 : : nodeCommMap,
1116 : : const std::map< int, std::vector< std::size_t > >& bface,
1117 : : const std::map< int, std::vector< std::size_t > >& bnode,
1118 : : const std::vector< std::size_t >& triinpoel )
1119 : : // *****************************************************************************
1120 : : // Receive new mesh from Refiner
1121 : : //! \param[in] ginpoel Mesh connectivity with global node ids
1122 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
1123 : : //! \param[in] coord New mesh node coordinates
1124 : : //! \param[in] addedNodes Newly added mesh nodes and their parents (local ids)
1125 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1126 : : //! \param[in] removedNodes Newly removed mesh node local ids
1127 : : //! \param[in] nodeCommMap New node communication map
1128 : : //! \param[in] bface Boundary-faces mapped to side set ids
1129 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
1130 : : //! \param[in] triinpoel Boundary-face connectivity
1131 : : // *****************************************************************************
1132 : : {
1133 : 0 : auto d = Disc();
1134 : :
1135 : 0 : d->Itf() = 0; // Zero field output iteration count if AMR
1136 : 0 : ++d->Itr(); // Increase number of iterations with a change in the mesh
1137 : :
1138 : : // Resize mesh data structures after mesh refinement
1139 : 0 : d->resizePostAMR( chunk, coord, nodeCommMap, removedNodes );
1140 : :
1141 : : Assert(coord[0].size() == m_u.nunk()-removedNodes.size()+addedNodes.size(),
1142 : : "Incorrect vector length post-AMR: expected length after resizing = " +
1143 : : std::to_string(coord[0].size()) + ", actual unknown vector length = " +
1144 : : std::to_string(m_u.nunk()-removedNodes.size()+addedNodes.size()));
1145 : :
1146 : : // Remove newly removed nodes from solution vectors
1147 : 0 : m_u.rm( removedNodes );
1148 : 0 : m_un.rm( removedNodes );
1149 : 0 : m_rhs.rm( removedNodes );
1150 : 0 : m_grad.rm( removedNodes );
1151 : :
1152 : : // Resize auxiliary solution vectors
1153 : : auto npoin = coord[0].size();
1154 : : m_u.resize( npoin );
1155 : : m_un.resize( npoin );
1156 : : m_rhs.resize( npoin );
1157 : : m_grad.resize( npoin );
1158 : :
1159 : : // Update solution on new mesh
1160 [ - - ]: 0 : for (const auto& n : addedNodes)
1161 [ - - ]: 0 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
1162 : : Assert(n.first < m_u.nunk(), "Added node index out of bounds post-AMR");
1163 : : Assert(n.second[0] < m_u.nunk() && n.second[1] < m_u.nunk(),
1164 : : "Indices of parent-edge nodes out of bounds post-AMR");
1165 : 0 : m_u(n.first,c) = (m_u(n.second[0],c) + m_u(n.second[1],c))/2.0;
1166 : : }
1167 : :
1168 : : // Update physical-boundary node-, face-, and element lists
1169 : : m_bnode = bnode;
1170 : : m_bface = bface;
1171 : 0 : m_triinpoel = tk::remap( triinpoel, d->Lid() );
1172 : :
1173 : 0 : auto meshid = d->MeshId();
1174 [ - - ]: 0 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1175 : 0 : CkCallback(CkReductionTarget(Transporter,resized), d->Tr()) );
1176 : 0 : }
1177 : :
1178 : : void
1179 : 1864 : RieCG::writeFields( CkCallback cb )
1180 : : // *****************************************************************************
1181 : : // Output mesh-based fields to file
1182 : : //! \param[in] cb Function to continue with after the write
1183 : : // *****************************************************************************
1184 : : {
1185 [ + + ]: 1864 : if (g_cfg.get< tag::benchmark >()) { cb.send(); return; }
1186 : :
1187 : 1280 : auto d = Disc();
1188 : : auto ncomp = m_u.nprop();
1189 : :
1190 : : // Field output
1191 : :
1192 : : std::vector< std::string > nodefieldnames
1193 [ - + ][ + + ]: 8960 : {"density", "xvelocity", "yvelocity", "zvelocity", "energy", "pressure"};
[ - + ][ - - ]
[ - - ]
1194 [ + - ]: 20 : if (g_cfg.get< tag::steady >()) nodefieldnames.push_back( "mach" );
1195 : :
1196 : : using tk::operator/=;
1197 [ + - ]: 1280 : auto r = m_u.extract(0);
1198 [ + - ][ + - ]: 1280 : auto u = m_u.extract(1); u /= r;
1199 [ + - ][ + - ]: 1280 : auto v = m_u.extract(2); v /= r;
1200 [ + - ][ + - ]: 1280 : auto w = m_u.extract(3); w /= r;
1201 [ + - ][ + - ]: 1280 : auto e = m_u.extract(4); e /= r;
1202 [ + - ][ + + ]: 1280 : std::vector< tk::real > pr( m_u.nunk() ), ma;
[ - - ]
1203 [ + + ][ + - ]: 1280 : if (g_cfg.get< tag::steady >()) ma.resize( m_u.nunk() );
1204 [ + + ]: 221497 : for (std::size_t i=0; i<pr.size(); ++i) {
1205 : 220217 : auto vv = u[i]*u[i] + v[i]*v[i] + w[i]*w[i];
1206 [ + + ]: 220217 : pr[i] = eos::pressure( r[i]*(e[i] - 0.5*vv) );
1207 [ + + ]: 220217 : if (g_cfg.get< tag::steady >()) {
1208 : 1114 : ma[i] = std::sqrt(vv) / eos::soundspeed( r[i], pr[i] );
1209 : : }
1210 : : }
1211 : :
1212 : : std::vector< std::vector< tk::real > > nodefields{
1213 : : std::move(r), std::move(u), std::move(v), std::move(w), std::move(e),
1214 [ - + ][ + + ]: 8960 : std::move(pr) };
[ + - ][ - - ]
[ - - ][ - - ]
1215 [ + + ]: 1280 : if (g_cfg.get< tag::steady >()) nodefields.push_back( std::move(ma) );
1216 : :
1217 [ + + ]: 1540 : for (std::size_t c=0; c<ncomp-5; ++c) {
1218 [ + - ]: 260 : nodefieldnames.push_back( "c" + std::to_string(c) );
1219 [ + - ]: 520 : nodefields.push_back( m_u.extract(5+c) );
1220 : : }
1221 : :
1222 : : // query function to evaluate analytic solution (if defined)
1223 [ + - ]: 1280 : auto sol = problems::SOL();
1224 : :
1225 [ + + ]: 1280 : if (sol) {
1226 : : const auto& coord = d->Coord();
1227 : : const auto& x = coord[0];
1228 : : const auto& y = coord[1];
1229 : : const auto& z = coord[2];
1230 : : auto an = m_u;
1231 [ + - ][ - - ]: 1103 : std::vector< tk::real > ap( m_u.nunk() );
1232 [ + + ]: 137879 : for (std::size_t i=0; i<an.nunk(); ++i) {
1233 [ - + ]: 273552 : auto s = sol( x[i], y[i], z[i], d->T() );
1234 : 136776 : s[1] /= s[0];
1235 : 136776 : s[2] /= s[0];
1236 : 136776 : s[3] /= s[0];
1237 : 136776 : s[4] /= s[0];
1238 [ + + ]: 909740 : for (std::size_t c=0; c<s.size(); ++c) an(i,c) = s[c];
1239 : 136776 : s[4] -= 0.5*(s[1]*s[1] + s[2]*s[2] + s[3]*s[3]);
1240 : 136776 : ap[i] = eos::pressure( s[0]*s[4] );
1241 : : }
1242 [ + + ]: 6618 : for (std::size_t c=0; c<5; ++c) {
1243 [ + - ]: 5515 : nodefieldnames.push_back( nodefieldnames[c] + "_analytic" );
1244 [ + - ]: 11030 : nodefields.push_back( an.extract(c) );
1245 : : }
1246 [ + - ][ + - ]: 2206 : nodefieldnames.push_back( nodefieldnames[5] + "_analytic" );
1247 : : nodefields.push_back( std::move(ap) );
1248 [ + + ]: 1339 : for (std::size_t c=0; c<ncomp-5; ++c) {
1249 [ + - ]: 236 : nodefieldnames.push_back( nodefieldnames[6+c] + "_analytic" );
1250 [ + - ][ - - ]: 472 : nodefields.push_back( an.extract(5+c) );
1251 : : }
1252 : : }
1253 : :
1254 : : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
1255 : :
1256 : : // Surface output
1257 : :
1258 : : std::vector< std::string > nodesurfnames;
1259 : : std::vector< std::vector< tk::real > > nodesurfs;
1260 : :
1261 : : const auto& f = g_cfg.get< tag::fieldout >();
1262 : :
1263 [ + + ]: 1280 : if (!f.empty()) {
1264 [ + - ]: 19 : nodesurfnames.push_back( "density" );
1265 [ + - ]: 19 : nodesurfnames.push_back( "xvelocity" );
1266 [ + - ]: 19 : nodesurfnames.push_back( "yvelocity" );
1267 [ + - ]: 19 : nodesurfnames.push_back( "zvelocity" );
1268 [ + - ]: 19 : nodesurfnames.push_back( "energy" );
1269 [ + - ]: 19 : nodesurfnames.push_back( "pressure" );
1270 : :
1271 [ - + ]: 19 : if (g_cfg.get< tag::steady >()) {
1272 [ - - ]: 0 : nodesurfnames.push_back( "mach" );
1273 : : }
1274 : :
1275 [ - + ]: 19 : for (std::size_t c=0; c<ncomp-5; ++c) {
1276 [ - - ]: 0 : nodesurfnames.push_back( "c" + std::to_string(c) );
1277 : : }
1278 : :
1279 [ + - ]: 19 : auto bnode = tk::bfacenodes( m_bface, m_triinpoel );
1280 [ + - ]: 19 : std::set< int > outsets( begin(f), end(f) );
1281 [ + + ]: 71 : for (auto sideset : outsets) {
1282 : : auto b = bnode.find(sideset);
1283 [ + + ]: 52 : if (b == end(bnode)) continue;
1284 : : const auto& nodes = b->second;
1285 : : auto i = nodesurfs.size();
1286 : 46 : auto ns = ncomp + 1;
1287 [ - + ]: 46 : if (g_cfg.get< tag::steady >()) ++ns;
1288 : : nodesurfs.insert( end(nodesurfs), ns,
1289 [ + - ]: 92 : std::vector< tk::real >( nodes.size() ) );
1290 : : std::size_t j = 0;
1291 [ + + ]: 4876 : for (auto n : nodes) {
1292 [ + - ]: 4830 : const auto s = m_u[n];
1293 : : std::size_t p = 0;
1294 : 4830 : nodesurfs[i+(p++)][j] = s[0];
1295 : 4830 : nodesurfs[i+(p++)][j] = s[1]/s[0];
1296 : 4830 : nodesurfs[i+(p++)][j] = s[2]/s[0];
1297 : 4830 : nodesurfs[i+(p++)][j] = s[3]/s[0];
1298 : 4830 : nodesurfs[i+(p++)][j] = s[4]/s[0];
1299 : 4830 : auto vv = (s[1]*s[1] + s[2]*s[2] + s[3]*s[3])/s[0]/s[0];
1300 : 4830 : auto ei = s[4]/s[0] - 0.5*vv;
1301 : 4830 : auto sp = eos::pressure( s[0]*ei );
1302 : 4830 : nodesurfs[i+(p++)][j] = sp;
1303 [ - + ]: 4830 : for (std::size_t c=0; c<ncomp-5; ++c) nodesurfs[i+(p++)+c][j] = s[5+c];
1304 [ - + ]: 4830 : if (g_cfg.get< tag::steady >()) {
1305 : 0 : nodesurfs[i+(p++)][j] = std::sqrt(vv) / eos::soundspeed( s[0], sp );
1306 : : }
1307 : 4830 : ++j;
1308 : : }
1309 : : }
1310 : : }
1311 : :
1312 : : // Send mesh and fields data (solution dump) for output to file
1313 [ + - ][ + - ]: 2560 : d->write( d->Inpoel(), d->Coord(), m_bface, tk::remap(m_bnode,d->Lid()),
1314 [ + - ]: 1280 : m_triinpoel, {}, nodefieldnames, {}, nodesurfnames,
1315 : : {}, nodefields, {}, nodesurfs, cb );
1316 [ + - ][ + - ]: 6400 : }
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ][ - - ]
[ - - ]
1317 : :
1318 : : void
1319 : 12680 : RieCG::out()
1320 : : // *****************************************************************************
1321 : : // Output mesh field data
1322 : : // *****************************************************************************
1323 : : {
1324 : 12680 : auto d = Disc();
1325 : :
1326 : : // Time history
1327 [ + + ][ + + ]: 12680 : if (d->histiter() or d->histtime() or d->histrange()) {
[ + + ]
1328 : : auto ncomp = m_u.nprop();
1329 : : const auto& inpoel = d->Inpoel();
1330 : 810 : std::vector< std::vector< tk::real > > hist( d->Hist().size() );
1331 : : std::size_t j = 0;
1332 [ + + ]: 1325 : for (const auto& p : d->Hist()) {
1333 [ + - ]: 515 : auto e = p.get< tag::elem >(); // host element id
1334 : : const auto& n = p.get< tag::fn >(); // shapefunctions evaluated at point
1335 [ + - ]: 515 : hist[j].resize( ncomp+1, 0.0 );
1336 [ + + ]: 2575 : for (std::size_t i=0; i<4; ++i) {
1337 [ + - ]: 2060 : const auto u = m_u[ inpoel[e*4+i] ];
1338 : 2060 : hist[j][0] += n[i] * u[0];
1339 : 2060 : hist[j][1] += n[i] * u[1]/u[0];
1340 : 2060 : hist[j][2] += n[i] * u[2]/u[0];
1341 : 2060 : hist[j][3] += n[i] * u[3]/u[0];
1342 : 2060 : hist[j][4] += n[i] * u[4]/u[0];
1343 : 2060 : auto ei = u[4]/u[0] - 0.5*(u[1]*u[1] + u[2]*u[2] + u[3]*u[3])/u[0]/u[0];
1344 : 2060 : hist[j][5] += n[i] * eos::pressure( u[0]*ei );
1345 [ - + ]: 2060 : for (std::size_t c=5; c<ncomp; ++c) hist[j][c+1] += n[i] * u[c];
1346 : : }
1347 : 515 : ++j;
1348 : : }
1349 [ + - ]: 810 : d->history( std::move(hist) );
1350 : 810 : }
1351 : :
1352 : : // Field data
1353 [ + + ][ + + ]: 12680 : if (d->fielditer() or d->fieldtime() or d->fieldrange() or m_finished) {
[ + + ][ + + ]
1354 [ + - ][ + - ]: 3915 : writeFields( CkCallback(CkIndex_RieCG::integrals(), thisProxy[thisIndex]) );
1355 : : } else {
1356 : 11375 : integrals();
1357 : : }
1358 : 12680 : }
1359 : :
1360 : : void
1361 : 12680 : RieCG::integrals()
1362 : : // *****************************************************************************
1363 : : // Compute integral quantities for output
1364 : : // *****************************************************************************
1365 : : {
1366 : 12680 : auto d = Disc();
1367 : :
1368 [ + + ][ + + ]: 12680 : if (d->integiter() or d->integtime() or d->integrange()) {
[ + + ]
1369 : :
1370 : : using namespace integrals;
1371 [ + - ]: 228 : std::vector< std::map< int, tk::real > > ints( NUMINT );
1372 : :
1373 : : // Prepend integral vector with metadata on the current time step:
1374 : : // current iteration count, current physical time, time step size
1375 [ + - ][ + - ]: 228 : ints[ ITER ][ 0 ] = static_cast< tk::real >( d->It() );
1376 [ + - ][ + - ]: 228 : ints[ TIME ][ 0 ] = d->T();
1377 [ + - ][ + - ]: 228 : ints[ DT ][ 0 ] = d->Dt();
1378 : :
1379 : : // Compute integrals requested for surfaces requested
1380 : : const auto& reqv = g_cfg.get< tag::integout_integrals >();
1381 : 228 : std::unordered_set< std::string > req( begin(reqv), end(reqv) );
1382 [ + - ][ + - ]: 456 : if (req.count("mass_flow_rate")) {
1383 [ + + ]: 376 : for (const auto& [s,sint] : m_surfint) {
1384 [ + - ]: 148 : auto& mfr = ints[ MASS_FLOW_RATE ][ s ];
1385 : : const auto& nodes = sint.first;
1386 : : const auto& ndA = sint.second;
1387 : : auto n = ndA.data();
1388 [ + + ]: 1776 : for (auto p : nodes) {
1389 : 1628 : mfr += n[0]*m_u(p,1) + n[1]*m_u(p,2) + n[2]*m_u(p,3);
1390 : 1628 : n += 3;
1391 : : }
1392 : : }
1393 : : }
1394 : :
1395 [ + - ]: 228 : auto stream = serialize( d->MeshId(), ints );
1396 [ + - ][ + - ]: 456 : d->contribute( stream.first, stream.second.get(), IntegralsMerger,
1397 [ + - ][ - - ]: 228 : CkCallback(CkIndex_Transporter::integrals(nullptr), d->Tr()) );
1398 : :
1399 : 228 : } else {
1400 : :
1401 : 12452 : step();
1402 : :
1403 : : }
1404 : 12680 : }
1405 : :
1406 : : void
1407 : 38040 : RieCG::stage()
1408 : : // *****************************************************************************
1409 : : // Evaluate whether to continue with next time step stage
1410 : : // *****************************************************************************
1411 : : {
1412 : : // Increment Runge-Kutta stage counter
1413 : 38040 : ++m_stage;
1414 : :
1415 : : // If not all Runge-Kutta stages complete, continue to next time stage,
1416 : : // otherwise output field data to file(s)
1417 [ + + ]: 38040 : if (m_stage < 3) grad(); else out();
1418 : 38040 : }
1419 : :
1420 : : void
1421 : 12121 : RieCG::evalLB( int nrestart )
1422 : : // *****************************************************************************
1423 : : // Evaluate whether to do load balancing
1424 : : //! \param[in] nrestart Number of times restarted
1425 : : // *****************************************************************************
1426 : : {
1427 : 12121 : auto d = Disc();
1428 : :
1429 : : // Detect if just returned from a checkpoint and if so, zero timers and
1430 : : // finished flag
1431 [ + + ]: 12121 : if (d->restarted( nrestart )) m_finished = 0;
1432 : :
1433 : 12121 : const auto lbfreq = g_cfg.get< tag::lbfreq >();
1434 [ + + ]: 12121 : const auto nonblocking = g_cfg.get< tag::nonblocking >();
1435 : :
1436 : : // Load balancing if user frequency is reached or after the second time-step
1437 [ + + ][ + + ]: 12121 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
1438 : :
1439 : 9039 : AtSync();
1440 [ - + ]: 9039 : if (nonblocking) dt();
1441 : :
1442 : : } else {
1443 : :
1444 : 3082 : dt();
1445 : :
1446 : : }
1447 : 12121 : }
1448 : :
1449 : : void
1450 : 12116 : RieCG::evalRestart()
1451 : : // *****************************************************************************
1452 : : // Evaluate whether to save checkpoint/restart
1453 : : // *****************************************************************************
1454 : : {
1455 : 12116 : auto d = Disc();
1456 : :
1457 : 12116 : const auto rsfreq = g_cfg.get< tag::rsfreq >();
1458 : 12116 : const auto benchmark = g_cfg.get< tag::benchmark >();
1459 : :
1460 [ + + ][ - + ]: 12116 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
1461 : :
1462 : 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
1463 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
1464 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
1465 : :
1466 : : } else {
1467 : :
1468 : 12116 : evalLB( /* nrestart = */ -1 );
1469 : :
1470 : : }
1471 : 12116 : }
1472 : :
1473 : : void
1474 : 12680 : RieCG::step()
1475 : : // *****************************************************************************
1476 : : // Evaluate whether to continue with next time step
1477 : : // *****************************************************************************
1478 : : {
1479 : 12680 : auto d = Disc();
1480 : :
1481 : : // Output one-liner status report to screen
1482 : 12680 : d->status();
1483 : : // Reset Runge-Kutta stage counter
1484 : 12680 : m_stage = 0;
1485 : :
1486 [ + + ]: 12680 : if (not m_finished) {
1487 : :
1488 : 12116 : evalRestart();
1489 : :
1490 : : } else {
1491 : :
1492 : 564 : auto meshid = d->MeshId();
1493 [ + - ]: 1128 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1494 : 564 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
1495 : :
1496 : : }
1497 : 12680 : }
1498 : :
1499 : : #include "NoWarning/riecg.def.h"
|
