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-2024 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 [ + - ]: 559 : 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 [ + - ]: 559 : m_bnode( bnode ),
60 [ + - ]: 559 : 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 [ + - ]: 559 : m_un( m_u.nunk(), m_u.nprop() ),
64 [ + - ]: 559 : 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 : 2236 : 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 : 559 : std::unordered_map< std::size_t, std::size_t > map;
84 : 559 : 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 [ + - ][ + + ]: 70658 : if (!map.count(p)) map[p] = n++;
[ + - ]
89 [ + + ]: 802602 : for (auto q : tk::Around(psup,p)) { // for each edge p-q
90 [ + - ][ + + ]: 731944 : if (!map.count(q)) map[q] = n++;
[ + - ]
91 : : }
92 : : }
93 : :
94 [ - + ][ - - ]: 559 : 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 : 559 : }
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 : 559 : std::unordered_map< int, std::unordered_set< std::size_t > > dir;
117 [ + + ]: 1553 : for (const auto& s : g_cfg.get< tag::bc_dir >()) {
118 [ + - ]: 994 : 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 [ + - ]: 994 : 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 [ + - ][ + - ]: 36994 : n.insert( tk::cref_find(lid,g) );
137 : : }
138 : : }
139 : : }
140 : :
141 : : // Collect unique set of nodes + Dirichlet BC components mask
142 : 559 : auto ncomp = m_u.nprop();
143 : 559 : auto nmask = ncomp + 1;
144 : 559 : const auto& dbc = g_cfg.get< tag::bc_dir >();
145 : 559 : 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 [ + - ]: 994 : auto n = dir.find( mask[0] );
149 [ + + ]: 994 : if (n != end(dir))
150 [ + + ]: 24739 : for (auto p : n->second) {
151 [ + - ]: 24142 : 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 : 559 : 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 [ + - ]: 32 : 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 : 559 : const auto& pbc_set = g_cfg.get< tag::bc_pre >();
183 [ + + ]: 559 : if (!pbc_set.empty()) {
184 : 16 : 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 : 16 : 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 : 559 : std::unordered_map< int, std::unordered_set< std::size_t > > sym;
209 [ + + ]: 1107 : for (auto s : g_cfg.get< tag::bc_sym >()) {
210 [ + - ]: 548 : 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 : 559 : std::unordered_map< int, std::unordered_set< std::size_t > > far;
223 [ + + ]: 563 : for (auto s : g_cfg.get< tag::bc_far >()) {
224 [ + - ]: 4 : 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 : 4874 : decltype(m_bnorm) exp;
268 [ + + ]: 36558 : for (auto i : nodes) {
269 [ + + ]: 103082 : for (const auto& [s,b] : m_bnorm) {
270 [ + - ]: 71398 : auto k = b.find(i);
271 [ + + ][ + - ]: 71398 : if (k != end(b)) exp[s][i] = k->second;
[ + - ]
272 : : }
273 : : }
274 [ + - ][ + - ]: 4874 : thisProxy[c].comnorm( exp );
275 : 4874 : }
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 : 559 : const auto& coord = d->Coord();
288 : 559 : const auto& gid = d->Gid();
289 : 559 : const auto& x = coord[0];
290 : 559 : const auto& y = coord[1];
291 : 559 : 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 : 67392 : auto n = tk::cross( ba, ca );
312 : 67392 : 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 : 134784 : (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 [ + - ]: 202176 : auto& v = m_bnorm[setid]; // associate side set id
324 [ + - ]: 202176 : 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 [ + - ]: 202176 : 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 : 559 : const auto& gid = d->Gid();
347 : 559 : const auto& inpoel = d->Inpoel();
348 : :
349 : 559 : const auto& coord = d->Coord();
350 : 559 : const auto& x = coord[0];
351 : 559 : const auto& y = coord[1];
352 : 559 : 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 : 1489254 : tk::real sig = 1.0;
371 [ + + ]: 1489254 : if (ed[0] > ed[1]) {
372 : 518294 : std::swap( ed[0], ed[1] );
373 : 518294 : sig = -1.0;
374 : : }
375 [ + - ]: 1489254 : 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 [ + - ]: 4257 : auto& bndnorm = m_bnormc[s];
393 [ + + ]: 19310 : for (const auto& [p,n] : b) {
394 [ + - ]: 15053 : 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 : 768 : 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 : 768 : NodeDiagnostics::registerReducers();
420 : 768 : IntegralsMerger = CkReduction::addReducer( integrals::mergeIntegrals );
421 : 768 : }
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 : 30 : auto ncomp = m_u.nprop();
457 [ + + ]: 34 : for (std::size_t c=5; c<ncomp; ++c)
458 [ + - ][ + - ]: 4 : 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 [ + - ]: 1354 : auto& bndnorm = m_bnorm[s];
493 [ + + ]: 13104 : for (const auto& [g,n] : b) {
494 [ + - ]: 11750 : 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 [ - + ][ - - ]: 47643 : 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 : 559 : decltype(m_bnorm) loc;
517 [ + + ]: 2089 : for (auto& [s,b] : m_bnorm) {
518 [ + - ]: 1530 : 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 : 559 : 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 : 559 : 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 : 559 : std::unordered_map< int, std::vector< std::size_t > > surfintnodes;
541 : 559 : 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 [ + - ]: 28 : 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 : 10 : nodes = std::move(n);
563 [ + - ]: 10 : ndA.resize( nodes.size()*3 );
564 : 10 : std::size_t a = 0;
565 [ + + ]: 120 : for (auto p : nodes) {
566 [ + - ]: 110 : 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 [ + - ]: 42000 : auto m = m_bnorm.find(s);
585 [ + + ]: 42000 : if (m != end(m_bnorm)) {
586 [ + - ]: 39485 : 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 [ + - ]: 48 : auto n = m_bnorm.find(s);
604 [ + + ]: 48 : if (n != end(m_bnorm)) {
605 [ + - ]: 24 : 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 [ - + ][ - - ]: 559 : Assert( !m_domedgeint.empty(), "No domain edges to group" );
[ - - ][ - - ]
628 : :
629 : : #ifndef NDEBUG
630 : 559 : 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 : 559 : tk::destroy( m_dsupedge[0] );
638 : 559 : tk::destroy( m_dsupedge[1] );
639 : 559 : tk::destroy( m_dsupedge[2] );
640 : :
641 : 559 : tk::destroy( m_dsupint[0] );
642 : 559 : tk::destroy( m_dsupint[1] );
643 : 559 : tk::destroy( m_dsupint[2] );
644 : :
645 : 559 : 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 : 248209 : int f = 0;
656 : : tk::real sig[6];
657 [ + - ][ + + ]: 1737463 : decltype(m_domedgeint)::const_iterator d[6];
658 [ + + ]: 711561 : for (const auto& [p,q] : tk::lpoed) {
659 : 677654 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
660 [ + + ]: 677654 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
661 [ + - ]: 677654 : d[f] = m_domedgeint.find( ed );
662 [ + + ]: 677654 : if (d[f] == end(m_domedgeint)) break; else ++f;
663 : : }
664 [ + + ]: 248209 : if (f == 6) {
665 [ + - ]: 33907 : m_dsupedge[0].push_back( N[0] );
666 [ + - ]: 33907 : m_dsupedge[0].push_back( N[1] );
667 [ + - ]: 33907 : m_dsupedge[0].push_back( N[2] );
668 [ + - ]: 33907 : m_dsupedge[0].push_back( N[3] );
669 [ + - ][ + + ]: 169535 : for (const auto& [a,b,c] : tk::lpofa) untri.erase( { N[a], N[b], N[c] } );
670 [ + + ]: 237349 : for (int ed=0; ed<6; ++ed) {
671 [ + - ]: 203442 : m_dsupint[0].push_back( sig[ed] * d[ed]->second[0] );
672 [ + - ]: 203442 : m_dsupint[0].push_back( sig[ed] * d[ed]->second[1] );
673 [ + - ]: 203442 : m_dsupint[0].push_back( sig[ed] * d[ed]->second[2] );
674 [ + - ]: 203442 : m_domedgeint.erase( d[ed] );
675 : : }
676 : : }
677 : : }
678 : :
679 [ + + ]: 409013 : for (const auto& N : untri) {
680 : 408454 : int f = 0;
681 : : tk::real sig[3];
682 [ + - ][ + + ]: 1633816 : decltype(m_domedgeint)::const_iterator d[3];
683 [ + + ]: 669696 : for (const auto& [p,q] : tk::lpoet) {
684 : 645171 : tk::UnsMesh::Edge ed{ gid[N[p]], gid[N[q]] };
685 [ + + ]: 645171 : sig[f] = ed[0] < ed[1] ? 1.0 : -1.0;
686 [ + - ]: 645171 : d[f] = m_domedgeint.find( ed );
687 [ + + ]: 645171 : if (d[f] == end(m_domedgeint)) break; else ++f;
688 : : }
689 [ + + ]: 408454 : if (f == 3) {
690 [ + - ]: 24525 : m_dsupedge[1].push_back( N[0] );
691 [ + - ]: 24525 : m_dsupedge[1].push_back( N[1] );
692 [ + - ]: 24525 : m_dsupedge[1].push_back( N[2] );
693 [ + + ]: 98100 : for (int ed=0; ed<3; ++ed) {
694 [ + - ]: 73575 : m_dsupint[1].push_back( sig[ed] * d[ed]->second[0] );
695 [ + - ]: 73575 : m_dsupint[1].push_back( sig[ed] * d[ed]->second[1] );
696 [ + - ]: 73575 : m_dsupint[1].push_back( sig[ed] * d[ed]->second[2] );
697 [ + - ]: 73575 : 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 : 559 : std::size_t k = 0;
705 [ + + ]: 89514 : for (const auto& [ed,d] : m_domedgeint) {
706 : 88955 : auto e = m_dsupedge[2].data() + k*2;
707 [ + - ]: 88955 : e[0] = tk::cref_find( lid, ed[0] );
708 [ + - ]: 88955 : e[1] = tk::cref_find( lid, ed[1] );
709 : 88955 : auto i = m_dsupint[2].data() + k*3;
710 : 88955 : i[0] = d[0];
711 : 88955 : i[1] = d[1];
712 : 88955 : i[2] = d[2];
713 : 88955 : ++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 [ - + ][ - - ]: 559 : 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 [ + - ][ + - ]: 559 : 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 : 12680 : 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 : 9760 : 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 [ + - ][ + - ]: 25360 : 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 : : // Set new time step size
860 [ + - ]: 12680 : if (m_stage == 0) Disc()->setdt( newdt );
861 : :
862 : 12680 : grad();
863 : 12680 : }
864 : :
865 : : void
866 : 38040 : RieCG::grad()
867 : : // *****************************************************************************
868 : : // Compute gradients for next time step
869 : : // *****************************************************************************
870 : : {
871 : 38040 : auto d = Disc();
872 : :
873 : 38040 : riemann::grad( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_u, m_grad );
874 : :
875 : : // Send gradient contributions to neighbor chares
876 [ + + ]: 38040 : if (d->NodeCommMap().empty()) {
877 : 1404 : comgrad_complete();
878 : : } else {
879 : 36636 : const auto& lid = d->Lid();
880 [ + + ]: 332910 : for (const auto& [c,n] : d->NodeCommMap()) {
881 : 296274 : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
882 [ + - ][ + - ]: 2119578 : for (auto g : n) exp[g] = m_grad[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
883 [ + - ][ + - ]: 296274 : thisProxy[c].comgrad( exp );
884 : 296274 : }
885 : : }
886 : 38040 : owngrad_complete();
887 : 38040 : }
888 : :
889 : : void
890 : 296274 : RieCG::comgrad(
891 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& ingrad )
892 : : // *****************************************************************************
893 : : // Receive contributions to node gradients on chare-boundaries
894 : : //! \param[in] ingrad Partial contributions to chare-boundary nodes. Key:
895 : : //! global mesh node IDs, value: contributions for all scalar components.
896 : : //! \details This function receives contributions to m_grad, which stores the
897 : : //! gradients at mesh nodes. While m_grad stores own contributions, m_gradc
898 : : //! collects the neighbor chare contributions during communication. This way
899 : : //! work on m_grad and m_gradc is overlapped. The two are combined in rhs().
900 : : // *****************************************************************************
901 : : {
902 : : using tk::operator+=;
903 [ + - ][ + - ]: 2119578 : for (const auto& [g,r] : ingrad) m_gradc[g] += r;
[ + + ]
904 : :
905 : : // When we have heard from all chares we communicate with, this chare is done
906 [ + + ]: 296274 : if (++m_ngrad == Disc()->NodeCommMap().size()) {
907 : 36636 : m_ngrad = 0;
908 : 36636 : comgrad_complete();
909 : : }
910 : 296274 : }
911 : :
912 : : void
913 : 38040 : RieCG::rhs()
914 : : // *****************************************************************************
915 : : // Compute right-hand side of transport equations
916 : : // *****************************************************************************
917 : : {
918 : 38040 : auto d = Disc();
919 : 38040 : const auto& lid = d->Lid();
920 : 38040 : const auto steady = g_cfg.get< tag::steady >();
921 : :
922 : : // Combine own and communicated contributions to gradients
923 [ + + ]: 1127304 : for (const auto& [g,r] : m_gradc) {
924 [ + - ]: 1089264 : auto i = tk::cref_find( lid, g );
925 [ + - ][ + + ]: 17931342 : for (std::size_t c=0; c<r.size(); ++c) m_grad(i,c) += r[c];
926 : : }
927 : 38040 : tk::destroy(m_gradc);
928 : :
929 : : // divide weak result in gradients by nodal volume
930 : 38040 : const auto& vol = d->Vol();
931 [ + + ]: 3614505 : for (std::size_t p=0; p<m_grad.nunk(); ++p)
932 [ + + ]: 60027930 : for (std::size_t c=0; c<m_grad.nprop(); ++c)
933 : 56451465 : m_grad(p,c) /= vol[p];
934 : :
935 : : // Compute own portion of right-hand side for all equations
936 [ + + ]: 38040 : auto prev_rkcoef = m_stage == 0 ? 0.0 : rkcoef[m_stage-1];
937 : :
938 [ + + ]: 38040 : if (steady) {
939 [ + + ]: 16860 : for (std::size_t p=0; p<m_tp.size(); ++p) m_tp[p] += prev_rkcoef * m_dtp[p];
940 : : }
941 : :
942 : 76080 : riemann::rhs( m_dsupedge, m_dsupint, d->Coord(), m_triinpoel, m_besym, m_grad,
943 : 76080 : m_u, d->V(), d->T(), m_tp, m_rhs );
944 : :
945 [ + + ]: 38040 : if (steady) {
946 [ + + ]: 16860 : for (std::size_t p=0; p<m_tp.size(); ++p) m_tp[p] -= prev_rkcoef * m_dtp[p];
947 : : }
948 : :
949 : : // Communicate rhs to other chares on chare-boundary
950 [ + + ]: 38040 : if (d->NodeCommMap().empty()) {
951 : 1404 : comrhs_complete();
952 : : } else {
953 [ + + ]: 332910 : for (const auto& [c,n] : d->NodeCommMap()) {
954 : 296274 : std::unordered_map< std::size_t, std::vector< tk::real > > exp;
955 [ + - ][ + - ]: 2119578 : for (auto g : n) exp[g] = m_rhs[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
956 [ + - ][ + - ]: 296274 : thisProxy[c].comrhs( exp );
957 : 296274 : }
958 : : }
959 : 38040 : ownrhs_complete();
960 : 38040 : }
961 : :
962 : : void
963 : 296274 : RieCG::comrhs(
964 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& inrhs )
965 : : // *****************************************************************************
966 : : // Receive contributions to right-hand side vector on chare-boundaries
967 : : //! \param[in] inrhs Partial contributions of RHS to chare-boundary nodes. Key:
968 : : //! global mesh node IDs, value: contributions for all scalar components.
969 : : //! \details This function receives contributions to m_rhs, which stores the
970 : : //! right hand side vector at mesh nodes. While m_rhs stores own
971 : : //! contributions, m_rhsc collects the neighbor chare contributions during
972 : : //! communication. This way work on m_rhs and m_rhsc is overlapped. The two
973 : : //! are combined in solve().
974 : : // *****************************************************************************
975 : : {
976 : : using tk::operator+=;
977 [ + - ][ + - ]: 2119578 : for (const auto& [g,r] : inrhs) m_rhsc[g] += r;
[ + + ]
978 : :
979 : : // When we have heard from all chares we communicate with, this chare is done
980 [ + + ]: 296274 : if (++m_nrhs == Disc()->NodeCommMap().size()) {
981 : 36636 : m_nrhs = 0;
982 : 36636 : comrhs_complete();
983 : : }
984 : 296274 : }
985 : :
986 : : void
987 : : // cppcheck-suppress unusedFunction
988 : 38040 : RieCG::solve()
989 : : // *****************************************************************************
990 : : // Advance systems of equations
991 : : // *****************************************************************************
992 : : {
993 [ + - ]: 38040 : auto d = Disc();
994 [ + - ]: 38040 : const auto lid = d->Lid();
995 : 38040 : const auto steady = g_cfg.get< tag::steady >();
996 : :
997 : : // Combine own and communicated contributions to rhs
998 [ + + ]: 1127304 : for (const auto& [g,r] : m_rhsc) {
999 [ + - ]: 1089264 : auto i = tk::cref_find( lid, g );
1000 [ + - ][ + + ]: 6703290 : for (std::size_t c=0; c<r.size(); ++c) m_rhs(i,c) += r[c];
1001 : : }
1002 : 38040 : tk::destroy(m_rhsc);
1003 : :
1004 : : // Update state at time n
1005 [ + + ][ + - ]: 38040 : if (m_stage == 0) m_un = m_u;
1006 : :
1007 : : // Advance solution
1008 : 38040 : auto dt = d->Dt();
1009 : 38040 : const auto& vol = d->Vol();
1010 [ + + ]: 3614505 : for (std::size_t i=0; i<m_u.nunk(); ++i) {
1011 [ + + ]: 3576465 : if (steady) dt = m_dtp[i];
1012 [ + + ]: 22393620 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
1013 [ + - ][ + - ]: 18817155 : m_u(i,c) = m_un(i,c) - rkcoef[m_stage] * dt * m_rhs(i,c) / vol[i];
[ + - ]
1014 : : }
1015 : : }
1016 : :
1017 : : // Configure and apply scalar source to solution (if defined)
1018 [ + - ]: 38040 : auto src = problems::PHYS_SRC();
1019 [ + + ][ + - ]: 38040 : if (src) src( d->Coord(), d->T(), m_u );
1020 : :
1021 : : // Enforce boundary conditions
1022 [ + - ]: 38040 : BC( d->T() + rkcoef[m_stage] * d->Dt() );
1023 : :
1024 [ + + ]: 38040 : if (m_stage < 2) {
1025 : :
1026 : : // Activate SDAG wait for next time step stage
1027 [ + - ][ + - ]: 25360 : thisProxy[ thisIndex ].wait4grad();
1028 [ + - ][ + - ]: 25360 : thisProxy[ thisIndex ].wait4rhs();
1029 : :
1030 : : // start next time step stage
1031 [ + - ]: 25360 : stage();
1032 : :
1033 : : } else {
1034 : :
1035 : : // Activate SDAG waits for finishing this time step stage
1036 [ + - ][ + - ]: 12680 : thisProxy[ thisIndex ].wait4stage();
1037 : : // Compute diagnostics, e.g., residuals
1038 : 12680 : auto diag_iter = g_cfg.get< tag::diag_iter >();
1039 [ + - ]: 12680 : auto diag = m_diag.rhocompute( *d, m_u, m_un, diag_iter );
1040 : : // Increase number of iterations and physical time
1041 [ + - ]: 12680 : d->next();
1042 : : // Advance physical time for local time stepping
1043 [ + + ]: 12680 : if (steady) {
1044 : : using tk::operator+=;
1045 [ + - ]: 50 : m_tp += m_dtp;
1046 : : }
1047 : : // Evaluate residuals
1048 [ + + ][ + - ]: 12680 : if (!diag) evalres( std::vector< tk::real >( m_u.nprop(), 1.0 ) );
[ + - ]
1049 : :
1050 : : }
1051 : 38040 : }
1052 : :
1053 : : void
1054 : 12680 : RieCG::evalres( const std::vector< tk::real >& l2res )
1055 : : // *****************************************************************************
1056 : : // Evaluate residuals
1057 : : //! \param[in] l2res L2-norms of the residual for each scalar component
1058 : : //! computed across the whole problem
1059 : : // *****************************************************************************
1060 : : {
1061 [ + + ]: 12680 : if (g_cfg.get< tag::steady >()) {
1062 : 50 : const auto rc = g_cfg.get< tag::rescomp >() - 1;
1063 : 50 : Disc()->residual( l2res[rc] );
1064 : : }
1065 : :
1066 : 12680 : refine();
1067 : 12680 : }
1068 : :
1069 : : void
1070 : 12680 : RieCG::refine()
1071 : : // *****************************************************************************
1072 : : // Optionally refine/derefine mesh
1073 : : // *****************************************************************************
1074 : : {
1075 : 12680 : auto d = Disc();
1076 : :
1077 : : // See if this is the last time step
1078 [ + + ]: 12680 : if (d->finished()) m_finished = 1;
1079 : :
1080 : 12680 : auto dtref = g_cfg.get< tag::href_dt >();
1081 : 12680 : auto dtfreq = g_cfg.get< tag::href_dtfreq >();
1082 : :
1083 : : // if t>0 refinement enabled and we hit the frequency
1084 [ - + ][ - - ]: 12680 : if (dtref && !(d->It() % dtfreq)) { // refine
[ - + ]
1085 : :
1086 : 0 : d->refined() = 1;
1087 : 0 : d->startvol();
1088 : 0 : d->Ref()->dtref( m_bface, m_bnode, m_triinpoel );
1089 : :
1090 : : // Activate SDAG waits for re-computing the integrals
1091 [ - - ][ - - ]: 0 : thisProxy[ thisIndex ].wait4int();
1092 : :
1093 : : } else { // do not refine
1094 : :
1095 : 12680 : d->refined() = 0;
1096 : 12680 : feop_complete();
1097 : 12680 : resize_complete();
1098 : :
1099 : : }
1100 : 12680 : }
1101 : :
1102 : : void
1103 : 0 : RieCG::resizePostAMR(
1104 : : const std::vector< std::size_t >& /*ginpoel*/,
1105 : : const tk::UnsMesh::Chunk& chunk,
1106 : : const tk::UnsMesh::Coords& coord,
1107 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
1108 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
1109 : : const std::set< std::size_t >& removedNodes,
1110 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
1111 : : nodeCommMap,
1112 : : const std::map< int, std::vector< std::size_t > >& bface,
1113 : : const std::map< int, std::vector< std::size_t > >& bnode,
1114 : : const std::vector< std::size_t >& triinpoel )
1115 : : // *****************************************************************************
1116 : : // Receive new mesh from Refiner
1117 : : //! \param[in] ginpoel Mesh connectivity with global node ids
1118 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
1119 : : //! \param[in] coord New mesh node coordinates
1120 : : //! \param[in] addedNodes Newly added mesh nodes and their parents (local ids)
1121 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1122 : : //! \param[in] removedNodes Newly removed mesh node local ids
1123 : : //! \param[in] nodeCommMap New node communication map
1124 : : //! \param[in] bface Boundary-faces mapped to side set ids
1125 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
1126 : : //! \param[in] triinpoel Boundary-face connectivity
1127 : : // *****************************************************************************
1128 : : {
1129 [ - - ]: 0 : auto d = Disc();
1130 : :
1131 : 0 : d->Itf() = 0; // Zero field output iteration count if AMR
1132 : 0 : ++d->Itr(); // Increase number of iterations with a change in the mesh
1133 : :
1134 : : // Resize mesh data structures after mesh refinement
1135 [ - - ]: 0 : d->resizePostAMR( chunk, coord, nodeCommMap, removedNodes );
1136 : :
1137 [ - - ][ - - ]: 0 : Assert(coord[0].size() == m_u.nunk()-removedNodes.size()+addedNodes.size(),
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
1138 : : "Incorrect vector length post-AMR: expected length after resizing = " +
1139 : : std::to_string(coord[0].size()) + ", actual unknown vector length = " +
1140 : : std::to_string(m_u.nunk()-removedNodes.size()+addedNodes.size()));
1141 : :
1142 : : // Remove newly removed nodes from solution vectors
1143 [ - - ]: 0 : m_u.rm( removedNodes );
1144 [ - - ]: 0 : m_un.rm( removedNodes );
1145 [ - - ]: 0 : m_rhs.rm( removedNodes );
1146 [ - - ]: 0 : m_grad.rm( removedNodes );
1147 : :
1148 : : // Resize auxiliary solution vectors
1149 : 0 : auto npoin = coord[0].size();
1150 [ - - ]: 0 : m_u.resize( npoin );
1151 [ - - ]: 0 : m_un.resize( npoin );
1152 [ - - ]: 0 : m_rhs.resize( npoin );
1153 [ - - ]: 0 : m_grad.resize( npoin );
1154 : :
1155 : : // Update solution on new mesh
1156 [ - - ]: 0 : for (const auto& n : addedNodes)
1157 [ - - ]: 0 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
1158 [ - - ][ - - ]: 0 : Assert(n.first < m_u.nunk(), "Added node index out of bounds post-AMR");
[ - - ][ - - ]
1159 [ - - ][ - - ]: 0 : Assert(n.second[0] < m_u.nunk() && n.second[1] < m_u.nunk(),
[ - - ][ - - ]
[ - - ]
1160 : : "Indices of parent-edge nodes out of bounds post-AMR");
1161 [ - - ][ - - ]: 0 : m_u(n.first,c) = (m_u(n.second[0],c) + m_u(n.second[1],c))/2.0;
[ - - ]
1162 : : }
1163 : :
1164 : : // Update physical-boundary node-, face-, and element lists
1165 [ - - ]: 0 : m_bnode = bnode;
1166 [ - - ]: 0 : m_bface = bface;
1167 [ - - ]: 0 : m_triinpoel = tk::remap( triinpoel, d->Lid() );
1168 : :
1169 : 0 : auto meshid = d->MeshId();
1170 [ - - ]: 0 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1171 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,resized), d->Tr()) );
1172 : 0 : }
1173 : :
1174 : : void
1175 : 1864 : RieCG::writeFields( CkCallback cb )
1176 : : // *****************************************************************************
1177 : : // Output mesh-based fields to file
1178 : : //! \param[in] cb Function to continue with after the write
1179 : : // *****************************************************************************
1180 : : {
1181 [ + + ][ + - ]: 1864 : if (g_cfg.get< tag::benchmark >()) { cb.send(); return; }
1182 : :
1183 [ + - ]: 1280 : auto d = Disc();
1184 : 1280 : auto ncomp = m_u.nprop();
1185 : :
1186 : : // Field output
1187 : :
1188 : : std::vector< std::string > nodefieldnames
1189 [ + - ][ + + ]: 8960 : {"density", "xvelocity", "yvelocity", "zvelocity", "energy", "pressure"};
[ - - ]
1190 [ + + ][ + - ]: 1280 : if (g_cfg.get< tag::steady >()) nodefieldnames.push_back( "mach" );
[ + - ]
1191 : :
1192 : : using tk::operator/=;
1193 [ + - ]: 1280 : auto r = m_u.extract(0);
1194 [ + - ][ + - ]: 1280 : auto u = m_u.extract(1); u /= r;
1195 [ + - ][ + - ]: 1280 : auto v = m_u.extract(2); v /= r;
1196 [ + - ][ + - ]: 1280 : auto w = m_u.extract(3); w /= r;
1197 [ + - ][ + - ]: 1280 : auto e = m_u.extract(4); e /= r;
1198 [ + - ]: 1280 : std::vector< tk::real > pr( m_u.nunk() ), ma;
1199 [ + + ][ + - ]: 1280 : if (g_cfg.get< tag::steady >()) ma.resize( m_u.nunk() );
1200 [ + + ]: 221497 : for (std::size_t i=0; i<pr.size(); ++i) {
1201 : 220217 : auto vv = u[i]*u[i] + v[i]*v[i] + w[i]*w[i];
1202 : 220217 : pr[i] = eos::pressure( r[i]*(e[i] - 0.5*vv) );
1203 [ + + ]: 220217 : if (g_cfg.get< tag::steady >()) {
1204 : 1114 : ma[i] = std::sqrt(vv) / eos::soundspeed( r[i], pr[i] );
1205 : : }
1206 : : }
1207 : :
1208 : : std::vector< std::vector< tk::real > > nodefields{
1209 : 6400 : std::move(r), std::move(u), std::move(v), std::move(w), std::move(e),
1210 [ + - ][ + + ]: 8960 : std::move(pr) };
[ - - ]
1211 [ + + ][ + - ]: 1280 : if (g_cfg.get< tag::steady >()) nodefields.push_back( std::move(ma) );
1212 : :
1213 [ + + ]: 1540 : for (std::size_t c=0; c<ncomp-5; ++c) {
1214 [ + - ][ + - ]: 260 : nodefieldnames.push_back( "c" + std::to_string(c) );
[ + - ]
1215 [ + - ][ + - ]: 260 : nodefields.push_back( m_u.extract(5+c) );
1216 : : }
1217 : :
1218 : : // query function to evaluate analytic solution (if defined)
1219 [ + - ]: 1280 : auto sol = problems::SOL();
1220 : :
1221 [ + + ]: 1280 : if (sol) {
1222 : 1103 : const auto& coord = d->Coord();
1223 : 1103 : const auto& x = coord[0];
1224 : 1103 : const auto& y = coord[1];
1225 : 1103 : const auto& z = coord[2];
1226 [ + - ]: 1103 : auto an = m_u;
1227 [ + - ]: 1103 : std::vector< tk::real > ap( m_u.nunk() );
1228 [ + + ]: 137879 : for (std::size_t i=0; i<an.nunk(); ++i) {
1229 [ + - ]: 136776 : auto s = sol( x[i], y[i], z[i], d->T() );
1230 : 136776 : s[1] /= s[0];
1231 : 136776 : s[2] /= s[0];
1232 : 136776 : s[3] /= s[0];
1233 : 136776 : s[4] /= s[0];
1234 [ + - ][ + + ]: 909740 : for (std::size_t c=0; c<s.size(); ++c) an(i,c) = s[c];
1235 : 136776 : s[4] -= 0.5*(s[1]*s[1] + s[2]*s[2] + s[3]*s[3]);
1236 : 136776 : ap[i] = eos::pressure( s[0]*s[4] );
1237 : 136776 : }
1238 [ + + ]: 6618 : for (std::size_t c=0; c<5; ++c) {
1239 [ + - ][ + - ]: 5515 : nodefieldnames.push_back( nodefieldnames[c] + "_analytic" );
1240 [ + - ][ + - ]: 5515 : nodefields.push_back( an.extract(c) );
1241 : : }
1242 [ + - ][ + - ]: 1103 : nodefieldnames.push_back( nodefieldnames[5] + "_analytic" );
1243 [ + - ]: 1103 : nodefields.push_back( std::move(ap) );
1244 [ + + ]: 1339 : for (std::size_t c=0; c<ncomp-5; ++c) {
1245 [ + - ][ + - ]: 236 : nodefieldnames.push_back( nodefieldnames[6+c] + "_analytic" );
1246 [ + - ][ + - ]: 236 : nodefields.push_back( an.extract(5+c) );
1247 : : }
1248 : 1103 : }
1249 : :
1250 [ - + ][ - - ]: 1280 : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
[ - - ][ - - ]
1251 : :
1252 : : // Surface output
1253 : :
1254 : 1280 : std::vector< std::string > nodesurfnames;
1255 : 1280 : std::vector< std::vector< tk::real > > nodesurfs;
1256 : :
1257 : 1280 : const auto& f = g_cfg.get< tag::fieldout >();
1258 : :
1259 [ + + ]: 1280 : if (!f.empty()) {
1260 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "density" );
1261 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "xvelocity" );
1262 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "yvelocity" );
1263 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "zvelocity" );
1264 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "energy" );
1265 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "pressure" );
1266 : :
1267 [ - + ]: 19 : if (g_cfg.get< tag::steady >()) {
1268 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "mach" );
1269 : : }
1270 : :
1271 [ - + ]: 19 : for (std::size_t c=0; c<ncomp-5; ++c) {
1272 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "c" + std::to_string(c) );
[ - - ]
1273 : : }
1274 : :
1275 [ + - ]: 19 : auto bnode = tk::bfacenodes( m_bface, m_triinpoel );
1276 [ + - ]: 19 : std::set< int > outsets( begin(f), end(f) );
1277 [ + + ]: 71 : for (auto sideset : outsets) {
1278 [ + - ]: 52 : auto b = bnode.find(sideset);
1279 [ + + ]: 52 : if (b == end(bnode)) continue;
1280 : 46 : const auto& nodes = b->second;
1281 : 46 : auto i = nodesurfs.size();
1282 : 46 : auto ns = ncomp + 1;
1283 [ - + ]: 46 : if (g_cfg.get< tag::steady >()) ++ns;
1284 [ + - ]: 46 : nodesurfs.insert( end(nodesurfs), ns,
1285 [ + - ]: 92 : std::vector< tk::real >( nodes.size() ) );
1286 : 46 : std::size_t j = 0;
1287 [ + + ]: 4876 : for (auto n : nodes) {
1288 [ + - ]: 4830 : const auto s = m_u[n];
1289 : 4830 : std::size_t p = 0;
1290 : 4830 : nodesurfs[i+(p++)][j] = s[0];
1291 : 4830 : nodesurfs[i+(p++)][j] = s[1]/s[0];
1292 : 4830 : nodesurfs[i+(p++)][j] = s[2]/s[0];
1293 : 4830 : nodesurfs[i+(p++)][j] = s[3]/s[0];
1294 : 4830 : nodesurfs[i+(p++)][j] = s[4]/s[0];
1295 : 4830 : auto vv = (s[1]*s[1] + s[2]*s[2] + s[3]*s[3])/s[0]/s[0];
1296 : 4830 : auto ei = s[4]/s[0] - 0.5*vv;
1297 : 4830 : auto sp = eos::pressure( s[0]*ei );
1298 : 4830 : nodesurfs[i+(p++)][j] = sp;
1299 [ - + ]: 4830 : for (std::size_t c=0; c<ncomp-5; ++c) nodesurfs[i+(p++)+c][j] = s[5+c];
1300 [ - + ]: 4830 : if (g_cfg.get< tag::steady >()) {
1301 : 0 : nodesurfs[i+(p++)][j] = std::sqrt(vv) / eos::soundspeed( s[0], sp );
1302 : : }
1303 : 4830 : ++j;
1304 : 4830 : }
1305 : : }
1306 : 19 : }
1307 : :
1308 : : // Send mesh and fields data (solution dump) for output to file
1309 [ + - ][ + - ]: 2560 : d->write( d->Inpoel(), d->Coord(), m_bface, tk::remap(m_bnode,d->Lid()),
1310 : 1280 : m_triinpoel, {}, nodefieldnames, {}, nodesurfnames,
1311 : : {}, nodefields, {}, nodesurfs, cb );
1312 [ + - ][ + - ]: 5120 : }
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
[ - - ][ - - ]
1313 : :
1314 : : void
1315 : 12680 : RieCG::out()
1316 : : // *****************************************************************************
1317 : : // Output mesh field data
1318 : : // *****************************************************************************
1319 : : {
1320 : 12680 : auto d = Disc();
1321 : :
1322 : : // Time history
1323 [ + + ][ + + ]: 12680 : if (d->histiter() or d->histtime() or d->histrange()) {
[ + + ][ + + ]
1324 : 810 : auto ncomp = m_u.nprop();
1325 : 810 : const auto& inpoel = d->Inpoel();
1326 [ + - ]: 810 : std::vector< std::vector< tk::real > > hist( d->Hist().size() );
1327 : 810 : std::size_t j = 0;
1328 [ + + ]: 1325 : for (const auto& p : d->Hist()) {
1329 : 515 : auto e = p.get< tag::elem >(); // host element id
1330 : 515 : const auto& n = p.get< tag::fn >(); // shapefunctions evaluated at point
1331 [ + - ]: 515 : hist[j].resize( ncomp+1, 0.0 );
1332 [ + + ]: 2575 : for (std::size_t i=0; i<4; ++i) {
1333 [ + - ]: 2060 : const auto u = m_u[ inpoel[e*4+i] ];
1334 : 2060 : hist[j][0] += n[i] * u[0];
1335 : 2060 : hist[j][1] += n[i] * u[1]/u[0];
1336 : 2060 : hist[j][2] += n[i] * u[2]/u[0];
1337 : 2060 : hist[j][3] += n[i] * u[3]/u[0];
1338 : 2060 : hist[j][4] += n[i] * u[4]/u[0];
1339 : 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];
1340 : 2060 : hist[j][5] += n[i] * eos::pressure( u[0]*ei );
1341 [ - + ]: 2060 : for (std::size_t c=5; c<ncomp; ++c) hist[j][c+1] += n[i] * u[c];
1342 : 2060 : }
1343 : 515 : ++j;
1344 : : }
1345 [ + - ]: 810 : d->history( std::move(hist) );
1346 : 810 : }
1347 : :
1348 : : // Field data
1349 [ + + ][ + + ]: 12680 : if (d->fielditer() or d->fieldtime() or d->fieldrange() or m_finished) {
[ + + ][ + + ]
[ + + ]
1350 [ + - ][ + - ]: 1305 : writeFields( CkCallback(CkIndex_RieCG::integrals(), thisProxy[thisIndex]) );
[ + - ][ + - ]
1351 : : } else {
1352 : 11375 : integrals();
1353 : : }
1354 : 12680 : }
1355 : :
1356 : : void
1357 : 12680 : RieCG::integrals()
1358 : : // *****************************************************************************
1359 : : // Compute integral quantities for output
1360 : : // *****************************************************************************
1361 : : {
1362 : 12680 : auto d = Disc();
1363 : :
1364 [ + + ][ + + ]: 12680 : if (d->integiter() or d->integtime() or d->integrange()) {
[ + + ][ + + ]
1365 : :
1366 : : using namespace integrals;
1367 [ + - ]: 228 : std::vector< std::map< int, tk::real > > ints( NUMINT );
1368 : : // Prepend integral vector with metadata on the current time step:
1369 : : // current iteration count, current physical time, time step size
1370 [ + - ]: 228 : ints[ ITER ][ 0 ] = static_cast< tk::real >( d->It() );
1371 [ + - ]: 228 : ints[ TIME ][ 0 ] = d->T();
1372 [ + - ]: 228 : ints[ DT ][ 0 ] = d->Dt();
1373 : : // Compute mass flow rate for surfaces requested
1374 [ + + ]: 376 : for (const auto& [s,sint] : m_surfint) {
1375 : : // cppcheck-suppress unreadVariable
1376 [ + - ]: 148 : auto& mfr = ints[ MASS_FLOW_RATE ][ s ];
1377 : 148 : const auto& nodes = sint.first;
1378 : 148 : const auto& ndA = sint.second;
1379 [ + + ]: 1776 : for (std::size_t i=0; i<nodes.size(); ++i) {
1380 : 1628 : auto p = nodes[i];
1381 [ + - ]: 1628 : mfr += ndA[i*3+0] * m_u(p,1)
1382 [ + - ]: 1628 : + ndA[i*3+1] * m_u(p,2)
1383 [ + - ]: 1628 : + ndA[i*3+2] * m_u(p,3);
1384 : : }
1385 : : }
1386 [ + - ]: 228 : auto stream = serialize( d->MeshId(), ints );
1387 [ + - ]: 228 : d->contribute( stream.first, stream.second.get(), IntegralsMerger,
1388 [ + - ][ + - ]: 456 : CkCallback(CkIndex_Transporter::integrals(nullptr), d->Tr()) );
1389 : :
1390 : 228 : } else {
1391 : :
1392 : 12452 : step();
1393 : :
1394 : : }
1395 : 12680 : }
1396 : :
1397 : : void
1398 : 38040 : RieCG::stage()
1399 : : // *****************************************************************************
1400 : : // Evaluate whether to continue with next time step stage
1401 : : // *****************************************************************************
1402 : : {
1403 : : // Increment Runge-Kutta stage counter
1404 : 38040 : ++m_stage;
1405 : :
1406 : : // If not all Runge-Kutta stages complete, continue to next time stage,
1407 : : // otherwise output field data to file(s)
1408 [ + + ]: 38040 : if (m_stage < 3) grad(); else out();
1409 : 38040 : }
1410 : :
1411 : : void
1412 : 12121 : RieCG::evalLB( int nrestart )
1413 : : // *****************************************************************************
1414 : : // Evaluate whether to do load balancing
1415 : : //! \param[in] nrestart Number of times restarted
1416 : : // *****************************************************************************
1417 : : {
1418 : 12121 : auto d = Disc();
1419 : :
1420 : : // Detect if just returned from a checkpoint and if so, zero timers and
1421 : : // finished flag
1422 [ + + ]: 12121 : if (d->restarted( nrestart )) m_finished = 0;
1423 : :
1424 : 12121 : const auto lbfreq = g_cfg.get< tag::lbfreq >();
1425 : 12121 : const auto nonblocking = g_cfg.get< tag::nonblocking >();
1426 : :
1427 : : // Load balancing if user frequency is reached or after the second time-step
1428 [ + + ][ + + ]: 12121 : if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
[ + + ]
1429 : :
1430 : 9039 : AtSync();
1431 [ - + ]: 9039 : if (nonblocking) dt();
1432 : :
1433 : : } else {
1434 : :
1435 : 3082 : dt();
1436 : :
1437 : : }
1438 : 12121 : }
1439 : :
1440 : : void
1441 : 12116 : RieCG::evalRestart()
1442 : : // *****************************************************************************
1443 : : // Evaluate whether to save checkpoint/restart
1444 : : // *****************************************************************************
1445 : : {
1446 : 12116 : auto d = Disc();
1447 : :
1448 : 12116 : const auto rsfreq = g_cfg.get< tag::rsfreq >();
1449 : 12116 : const auto benchmark = g_cfg.get< tag::benchmark >();
1450 : :
1451 [ + + ][ - + ]: 12116 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
[ - + ]
1452 : :
1453 [ - - ]: 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
1454 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
1455 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
1456 : :
1457 : 0 : } else {
1458 : :
1459 : 12116 : evalLB( /* nrestart = */ -1 );
1460 : :
1461 : : }
1462 : 12116 : }
1463 : :
1464 : : void
1465 : 12680 : RieCG::step()
1466 : : // *****************************************************************************
1467 : : // Evaluate whether to continue with next time step
1468 : : // *****************************************************************************
1469 : : {
1470 : 12680 : auto d = Disc();
1471 : :
1472 : : // Output one-liner status report to screen
1473 : 12680 : d->status();
1474 : : // Reset Runge-Kutta stage counter
1475 : 12680 : m_stage = 0;
1476 : :
1477 [ + + ]: 12680 : if (not m_finished) {
1478 : :
1479 : 12116 : evalRestart();
1480 : :
1481 : : } else {
1482 : :
1483 : 564 : auto meshid = d->MeshId();
1484 [ + - ]: 564 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1485 [ + - ][ + - ]: 1128 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
1486 : :
1487 : : }
1488 : 12680 : }
1489 : :
1490 : : #include "NoWarning/riecg.def.h"
|
