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