Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Inciter/KozCG.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 KozCG: Taylor-Galerkin, FCT, element-based continuous Galerkin
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "XystBuildConfig.hpp"
14 : : #include "KozCG.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 "Kozak.hpp"
30 : : #include "Problems.hpp"
31 : : #include "EOS.hpp"
32 : : #include "BC.hpp"
33 : :
34 : : namespace inciter {
35 : :
36 : : extern ctr::Config g_cfg;
37 : :
38 : : static CkReduction::reducerType IntegralsMerger;
39 : :
40 : : } // inciter::
41 : :
42 : : using inciter::g_cfg;
43 : : using inciter::KozCG;
44 : :
45 : 459 : KozCG::KozCG( const CProxy_Discretization& disc,
46 : : const std::map< int, std::vector< std::size_t > >& bface,
47 : : const std::map< int, std::vector< std::size_t > >& bnode,
48 : 459 : const std::vector< std::size_t >& triinpoel ) :
49 [ + - ]: 459 : m_disc( disc ),
50 : 459 : m_nrhs( 0 ),
51 : 459 : m_nnorm( 0 ),
52 : 459 : m_naec( 0 ),
53 : 459 : m_nalw( 0 ),
54 : 459 : m_nlim( 0 ),
55 [ + - ]: 459 : m_bnode( bnode ),
56 [ + - ]: 459 : m_bface( bface ),
57 [ + - ][ + - ]: 459 : m_triinpoel( tk::remap( triinpoel, Disc()->Lid() ) ),
58 [ + - ][ + - ]: 459 : m_u( Disc()->Gid().size(), g_cfg.get< tag::problem_ncomp >() ),
59 [ + - ]: 459 : m_p( m_u.nunk(), m_u.nprop()*2 ),
60 [ + - ]: 459 : m_q( m_u.nunk(), m_u.nprop()*2 ),
61 [ + - ]: 459 : m_a( m_u.nunk(), m_u.nprop() ),
62 [ + - ]: 459 : m_rhs( m_u.nunk(), m_u.nprop() ),
63 [ + - ]: 459 : m_dtp( m_u.nunk(), 0.0 ),
64 [ + - ]: 459 : m_tp( m_u.nunk(), g_cfg.get< tag::t0 >() ),
65 : 459 : m_finished( 0 ),
66 : 2295 : m_freezeflow( 1.0 )
67 : : // *****************************************************************************
68 : : // Constructor
69 : : //! \param[in] disc Discretization proxy
70 : : //! \param[in] bface Boundary-faces mapped to side sets used in the input file
71 : : //! \param[in] bnode Boundary-node lists mapped to side sets used in input file
72 : : //! \param[in] triinpoel Boundary-face connectivity where BCs set (global ids)
73 : : // *****************************************************************************
74 : : {
75 : 459 : usesAtSync = true; // enable migration at AtSync
76 : :
77 [ + - ]: 459 : auto d = Disc();
78 : :
79 : : // Compute total box IC volume
80 [ + - ]: 459 : d->boxvol();
81 : :
82 : : // Activate SDAG wait for initially computing integrals
83 [ + - ][ + - ]: 459 : thisProxy[ thisIndex ].wait4int();
84 : 459 : }
85 : :
86 : : void
87 : 459 : KozCG::setupBC()
88 : : // *****************************************************************************
89 : : // Prepare boundary condition data structures
90 : : // *****************************************************************************
91 : : {
92 : : // Query Dirichlet BC nodes associated to side sets
93 : 459 : std::unordered_map< int, std::unordered_set< std::size_t > > dir;
94 [ + + ]: 1161 : for (const auto& s : g_cfg.get< tag::bc_dir >()) {
95 [ + - ]: 702 : auto k = m_bface.find(s[0]);
96 [ + + ]: 702 : if (k != end(m_bface)) {
97 [ + - ]: 371 : auto& n = dir[ k->first ];
98 [ + + ]: 14233 : for (auto f : k->second) {
99 : 13862 : const auto t = m_triinpoel.data() + f*3;
100 [ + - ]: 13862 : n.insert( t[0] );
101 [ + - ]: 13862 : n.insert( t[1] );
102 [ + - ]: 13862 : n.insert( t[2] );
103 : : }
104 : : }
105 : : }
106 : :
107 : : // Augment Dirichlet BC nodes with nodes not necessarily part of faces
108 [ + - ]: 459 : const auto& lid = Disc()->Lid();
109 [ + + ]: 1161 : for (const auto& s : g_cfg.get< tag::bc_dir >()) {
110 [ + - ]: 702 : auto k = m_bnode.find(s[0]);
111 [ + + ]: 702 : if (k != end(m_bnode)) {
112 [ + - ]: 385 : auto& n = dir[ k->first ];
113 [ + + ]: 22105 : for (auto g : k->second) {
114 [ + - ][ + - ]: 21720 : n.insert( tk::cref_find(lid,g) );
115 : : }
116 : : }
117 : : }
118 : :
119 : : // Collect unique set of nodes + Dirichlet BC components mask
120 : 459 : auto ncomp = m_u.nprop();
121 : 459 : auto nmask = ncomp + 1;
122 : 459 : const auto& dbc = g_cfg.get< tag::bc_dir >();
123 : 459 : std::unordered_map< std::size_t, std::vector< int > > dirbcset;
124 [ + + ]: 1161 : for (const auto& mask : dbc) {
125 [ - + ][ - - ]: 702 : ErrChk( mask.size() == nmask, "Incorrect Dirichlet BC mask ncomp" );
[ - - ][ - - ]
126 [ + - ]: 702 : auto n = dir.find( mask[0] );
127 [ + + ]: 702 : if (n != end(dir)) {
128 [ + + ]: 10910 : for (auto p : n->second) {
129 [ + - ]: 10525 : auto& m = dirbcset[p];
130 [ + + ][ + - ]: 10525 : if (m.empty()) m.resize( ncomp, 0 );
131 [ + + ]: 68245 : for (std::size_t c=0; c<ncomp; ++c) {
132 [ + + ]: 57720 : if (!m[c]) m[c] = mask[c+1]; // overwrite mask if 0 -> 1
133 : : }
134 : : }
135 : : }
136 : : }
137 : : // Compile streamable list of nodes + Dirichlet BC components mask
138 : 459 : tk::destroy( m_dirbcmasks );
139 [ + + ]: 9298 : for (const auto& [p,mask] : dirbcset) {
140 [ + - ]: 8839 : m_dirbcmasks.push_back( p );
141 [ + - ]: 8839 : m_dirbcmasks.insert( end(m_dirbcmasks), begin(mask), end(mask) );
142 : : }
143 [ - + ][ - - ]: 459 : ErrChk( m_dirbcmasks.size() % nmask == 0, "Dirichlet BC masks incomplete" );
[ - - ][ - - ]
144 : :
145 : : // Query pressure BC nodes associated to side sets
146 : 459 : std::unordered_map< int, std::unordered_set< std::size_t > > pre;
147 [ - + ]: 459 : for (const auto& ss : g_cfg.get< tag::bc_pre >()) {
148 [ - - ]: 0 : for (const auto& s : ss) {
149 [ - - ]: 0 : auto k = m_bface.find(s);
150 [ - - ]: 0 : if (k != end(m_bface)) {
151 [ - - ]: 0 : auto& n = pre[ k->first ];
152 [ - - ]: 0 : for (auto f : k->second) {
153 : 0 : const auto t = m_triinpoel.data() + f*3;
154 [ - - ]: 0 : n.insert( t[0] );
155 [ - - ]: 0 : n.insert( t[1] );
156 [ - - ]: 0 : n.insert( t[2] );
157 : : }
158 : : }
159 : : }
160 : : }
161 : :
162 : : // Augment Pressure BC nodes with nodes not necessarily part of faces
163 [ - + ]: 459 : for (const auto& s : g_cfg.get< tag::bc_pre >()) {
164 [ - - ]: 0 : auto k = m_bnode.find(s[0]);
165 [ - - ]: 0 : if (k != end(m_bnode)) {
166 [ - - ]: 0 : auto& n = pre[ k->first ];
167 [ - - ]: 0 : for (auto g : k->second) {
168 [ - - ][ - - ]: 0 : n.insert( tk::cref_find(lid,g) );
169 : : }
170 : : }
171 : : }
172 : :
173 : : // Prepare density and pressure values for pressure BC nodes
174 : 459 : const auto& pbc_set = g_cfg.get< tag::bc_pre >();
175 [ - + ]: 459 : if (!pbc_set.empty()) {
176 : 0 : const auto& pbc_r = g_cfg.get< tag::bc_pre_density >();
177 [ - - ][ - - ]: 0 : ErrChk( pbc_r.size() == pbc_set.size(), "Pressure BC density unspecified" );
[ - - ][ - - ]
178 : 0 : const auto& pbc_p = g_cfg.get< tag::bc_pre_pressure >();
179 [ - - ][ - - ]: 0 : ErrChk( pbc_p.size() == pbc_set.size(), "Pressure BC pressure unspecified" );
[ - - ][ - - ]
180 : 0 : tk::destroy( m_prebcnodes );
181 : 0 : tk::destroy( m_prebcvals );
182 [ - - ]: 0 : for (const auto& [s,n] : pre) {
183 [ - - ]: 0 : m_prebcnodes.insert( end(m_prebcnodes), begin(n), end(n) );
184 [ - - ]: 0 : for (std::size_t p=0; p<pbc_set.size(); ++p) {
185 [ - - ]: 0 : for (auto u : pbc_set[p]) {
186 [ - - ]: 0 : if (s == u) {
187 [ - - ]: 0 : for (std::size_t i=0; i<n.size(); ++i) {
188 [ - - ]: 0 : m_prebcvals.push_back( pbc_r[p] );
189 [ - - ]: 0 : m_prebcvals.push_back( pbc_p[p] );
190 : : }
191 : : }
192 : : }
193 : : }
194 : : }
195 [ - - ][ - - ]: 0 : ErrChk( m_prebcnodes.size()*2 == m_prebcvals.size(),
[ - - ][ - - ]
196 : : "Pressure BC data incomplete" );
197 : : }
198 : :
199 : : // Query symmetry BC nodes associated to side sets
200 : 459 : std::unordered_map< int, std::unordered_set< std::size_t > > sym;
201 [ + + ]: 855 : for (auto s : g_cfg.get< tag::bc_sym >()) {
202 [ + - ]: 396 : auto k = m_bface.find(s);
203 [ + + ]: 396 : if (k != end(m_bface)) {
204 [ + - ]: 214 : auto& n = sym[ k->first ];
205 [ + + ]: 6870 : for (auto f : k->second) {
206 : 6656 : const auto t = m_triinpoel.data() + f*3;
207 [ + - ]: 6656 : n.insert( t[0] );
208 [ + - ]: 6656 : n.insert( t[1] );
209 [ + - ]: 6656 : n.insert( t[2] );
210 : : }
211 : : }
212 : : }
213 : :
214 : : // Query farfield BC nodes associated to side sets
215 : 459 : std::unordered_map< int, std::unordered_set< std::size_t > > far;
216 [ - + ]: 459 : for (auto s : g_cfg.get< tag::bc_far >()) {
217 [ - - ]: 0 : auto k = m_bface.find(s);
218 [ - - ]: 0 : if (k != end(m_bface)) {
219 [ - - ]: 0 : auto& n = far[ k->first ];
220 [ - - ]: 0 : for (auto f : k->second) {
221 : 0 : const auto t = m_triinpoel.data() + f*3;
222 [ - - ]: 0 : n.insert( t[0] );
223 [ - - ]: 0 : n.insert( t[1] );
224 [ - - ]: 0 : n.insert( t[2] );
225 : : }
226 : : }
227 : : }
228 : :
229 : : // Generate unique set of symmetry BC nodes
230 : 459 : tk::destroy( m_symbcnodeset );
231 [ + - ][ + + ]: 673 : for (const auto& [s,n] : sym) m_symbcnodeset.insert( begin(n), end(n) );
232 : : // Generate unique set of farfield BC nodes
233 : 459 : tk::destroy( m_farbcnodeset );
234 [ - - ][ - + ]: 459 : for (const auto& [s,n] : far) m_farbcnodeset.insert( begin(n), end(n) );
235 : :
236 : : // If farfield BC is set on a node, will not also set symmetry BC
237 [ - - ][ - + ]: 459 : for (auto i : m_farbcnodeset) m_symbcnodeset.erase(i);
238 : : // If pressure BC is set on a node, will not also set symmetry BC
239 [ - - ][ - + ]: 459 : for (auto i : m_prebcnodes) m_symbcnodeset.erase(i);
240 : 459 : }
241 : :
242 : : void
243 : 459 : KozCG::feop()
244 : : // *****************************************************************************
245 : : // Start (re-)computing finite element domain and boundary operators
246 : : // *****************************************************************************
247 : : {
248 : 459 : auto d = Disc();
249 : :
250 : : // Prepare boundary conditions data structures
251 : 459 : setupBC();
252 : :
253 : : // Compute local contributions to boundary normals and integrals
254 : 459 : bndint();
255 : :
256 : : // Send boundary point normal contributions to neighbor chares
257 [ + + ]: 459 : if (d->NodeCommMap().empty()) {
258 : 14 : comnorm_complete();
259 : : } else {
260 [ + + ]: 4923 : for (const auto& [c,nodes] : d->NodeCommMap()) {
261 : 4478 : decltype(m_bnorm) exp;
262 [ + + ]: 24342 : for (auto i : nodes) {
263 [ + + ]: 61843 : for (const auto& [s,b] : m_bnorm) {
264 [ + - ]: 41979 : auto k = b.find(i);
265 [ + + ][ + - ]: 41979 : if (k != end(b)) exp[s][i] = k->second;
[ + - ]
266 : : }
267 : : }
268 [ + - ][ + - ]: 4478 : thisProxy[c].comnorm( exp );
269 : 4478 : }
270 : : }
271 : 459 : ownnorm_complete();
272 : 459 : }
273 : :
274 : : void
275 : 459 : KozCG::bndint()
276 : : // *****************************************************************************
277 : : //! Compute local contributions to boundary normals and integrals
278 : : // *****************************************************************************
279 : : {
280 [ + - ]: 459 : auto d = Disc();
281 : 459 : const auto& coord = d->Coord();
282 : 459 : const auto& gid = d->Gid();
283 : 459 : const auto& x = coord[0];
284 : 459 : const auto& y = coord[1];
285 : 459 : const auto& z = coord[2];
286 : :
287 : : // Lambda to compute the inverse distance squared between boundary face
288 : : // centroid and boundary point. Here p is the global node id and c is the
289 : : // the boundary face centroid.
290 : 79914 : auto invdistsq = [&]( const tk::real c[], std::size_t p ){
291 : 79914 : return 1.0 / ( (c[0] - x[p]) * (c[0] - x[p]) +
292 : 79914 : (c[1] - y[p]) * (c[1] - y[p]) +
293 : 79914 : (c[2] - z[p]) * (c[2] - z[p]) );
294 : 459 : };
295 : :
296 : 459 : tk::destroy( m_bnorm );
297 : 459 : tk::destroy( m_bndpoinint );
298 : :
299 [ + + ]: 1544 : for (const auto& [ setid, faceids ] : m_bface) { // for all side sets
300 [ + + ]: 27723 : for (auto f : faceids) { // for all side set triangles
301 : : const std::array< std::size_t, 3 >
302 : 26638 : N{ m_triinpoel[f*3+0], m_triinpoel[f*3+1], m_triinpoel[f*3+2] };
303 : : const std::array< tk::real, 3 >
304 : 26638 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
305 : 26638 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
306 : 26638 : auto n = tk::cross( ba, ca );
307 : 26638 : auto A = tk::length( n );
308 : 26638 : n[0] /= A;
309 : 26638 : n[1] /= A;
310 : 26638 : n[2] /= A;
311 : 26638 : A /= 2.0;
312 : : const tk::real centroid[3] = {
313 : 26638 : (x[N[0]] + x[N[1]] + x[N[2]]) / 3.0,
314 : 26638 : (y[N[0]] + y[N[1]] + y[N[2]]) / 3.0,
315 : 53276 : (z[N[0]] + z[N[1]] + z[N[2]]) / 3.0 };
316 : :
317 : : // contribute all edges of triangle
318 [ + + ]: 106552 : for (const auto& [i,j] : tk::lpoet) {
319 : 79914 : auto p = N[i];
320 : 79914 : tk::real r = invdistsq( centroid, p );
321 [ + - ]: 79914 : auto& v = m_bnorm[setid]; // associate side set id
322 [ + - ]: 79914 : auto& bn = v[gid[p]]; // associate global node id of bnd pnt
323 : 79914 : bn[0] += r * n[0]; // inv.dist.sq-weighted normal
324 : 79914 : bn[1] += r * n[1];
325 : 79914 : bn[2] += r * n[2];
326 : 79914 : bn[3] += r; // inv.dist.sq of node from centroid
327 [ + - ]: 79914 : auto& b = m_bndpoinint[gid[p]];// assoc global id of bnd point
328 : 79914 : b[0] += n[0] * A / 3.0; // bnd-point integral
329 : 79914 : b[1] += n[1] * A / 3.0;
330 : 79914 : b[2] += n[2] * A / 3.0;
331 : : }
332 : : }
333 : : }
334 : 459 : }
335 : :
336 : : void
337 : 4478 : KozCG::comnorm( const decltype(m_bnorm)& inbnd )
338 : : // *****************************************************************************
339 : : // Receive contributions to boundary point normals on chare-boundaries
340 : : //! \param[in] inbnd Incoming partial sums of boundary point normals
341 : : // *****************************************************************************
342 : : {
343 : : // Buffer up incoming boundary point normal vector contributions
344 [ + + ]: 8116 : for (const auto& [s,b] : inbnd) {
345 [ + - ]: 3638 : auto& bndnorm = m_bnormc[s];
346 [ + + ]: 14967 : for (const auto& [p,n] : b) {
347 [ + - ]: 11329 : auto& norm = bndnorm[p];
348 : 11329 : norm[0] += n[0];
349 : 11329 : norm[1] += n[1];
350 : 11329 : norm[2] += n[2];
351 : 11329 : norm[3] += n[3];
352 : : }
353 : : }
354 : :
355 [ + + ]: 4478 : if (++m_nnorm == Disc()->NodeCommMap().size()) {
356 : 445 : m_nnorm = 0;
357 : 445 : comnorm_complete();
358 : : }
359 : 4478 : }
360 : :
361 : : void
362 : 768 : KozCG::registerReducers()
363 : : // *****************************************************************************
364 : : // Configure Charm++ reduction types initiated from this chare array
365 : : //! \details Since this is a [initnode] routine, the runtime system executes the
366 : : //! routine exactly once on every logical node early on in the Charm++ init
367 : : //! sequence. Must be static as it is called without an object. See also:
368 : : //! Section "Initializations at Program Startup" at in the Charm++ manual
369 : : //! http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
370 : : // *****************************************************************************
371 : : {
372 : 768 : NodeDiagnostics::registerReducers();
373 : 768 : IntegralsMerger = CkReduction::addReducer( integrals::mergeIntegrals );
374 : 768 : }
375 : :
376 : : void
377 : : // cppcheck-suppress unusedFunction
378 : 5549 : KozCG::ResumeFromSync()
379 : : // *****************************************************************************
380 : : // Return from migration
381 : : //! \details This is called when load balancing (LB) completes. The presence of
382 : : //! this function does not affect whether or not we block on LB.
383 : : // *****************************************************************************
384 : : {
385 [ - + ][ - - ]: 5549 : if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
[ - - ][ - - ]
386 : :
387 [ + - ]: 5549 : if (!g_cfg.get< tag::nonblocking >()) dt();
388 : 5549 : }
389 : :
390 : : void
391 : 459 : KozCG::setup( tk::real v )
392 : : // *****************************************************************************
393 : : // Start setup for solution
394 : : //! \param[in] v Total volume within user-specified box
395 : : // *****************************************************************************
396 : : {
397 : 459 : auto d = Disc();
398 : :
399 : : // Store user-defined box IC volume
400 : 459 : d->Boxvol() = v;
401 : :
402 : : // Set initial conditions
403 : 459 : problems::initialize( d->Coord(), m_u, d->T(), d->BoxNodes() );
404 : :
405 : : // Query time history field output labels from all PDEs integrated
406 [ + + ]: 459 : if (!g_cfg.get< tag::histout >().empty()) {
407 : : std::vector< std::string > var
408 [ + - ][ + + ]: 154 : {"density", "xvelocity", "yvelocity", "zvelocity", "energy", "pressure"};
[ - - ]
409 : 22 : auto ncomp = m_u.nprop();
410 [ - + ]: 22 : for (std::size_t c=5; c<ncomp; ++c)
411 [ - - ][ - - ]: 0 : var.push_back( "c" + std::to_string(c-5) );
[ - - ]
412 [ + - ]: 22 : d->histheader( std::move(var) );
413 : 22 : }
414 : :
415 : : // Compute finite element operators
416 : 459 : feop();
417 [ + - ][ + - ]: 503 : }
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
418 : :
419 : : void
420 : 459 : KozCG::start()
421 : : // *****************************************************************************
422 : : // Start time stepping
423 : : // *****************************************************************************
424 : : {
425 : : // Set flag that indicates that we are now during time stepping
426 : 459 : Disc()->Initial( 0 );
427 : : // Start timer measuring time stepping wall clock time
428 : 459 : Disc()->Timer().zero();
429 : : // Zero grind-timer
430 : 459 : Disc()->grindZero();
431 : : // Continue to first time step
432 : 459 : dt();
433 : 459 : }
434 : :
435 : : void
436 : 459 : KozCG::bnorm()
437 : : // *****************************************************************************
438 : : // Combine own and communicated portions of the boundary point normals
439 : : // *****************************************************************************
440 : : {
441 [ + - ]: 459 : const auto& lid = Disc()->Lid();
442 : :
443 : : // Combine own and communicated contributions to boundary point normals
444 [ + + ]: 1522 : for (const auto& [s,b] : m_bnormc) {
445 [ + - ]: 1063 : auto& bndnorm = m_bnorm[s];
446 [ + + ]: 9423 : for (const auto& [g,n] : b) {
447 [ + - ]: 8360 : auto& norm = bndnorm[g];
448 : 8360 : norm[0] += n[0];
449 : 8360 : norm[1] += n[1];
450 : 8360 : norm[2] += n[2];
451 : 8360 : norm[3] += n[3];
452 : : }
453 : : }
454 : 459 : tk::destroy( m_bnormc );
455 : :
456 : : // Divide summed point normals by the sum of the inverse distance squared
457 [ + + ]: 1612 : for (auto& [s,b] : m_bnorm) {
458 [ + + ]: 23016 : for (auto& [g,n] : b) {
459 : 21863 : n[0] /= n[3];
460 : 21863 : n[1] /= n[3];
461 : 21863 : n[2] /= n[3];
462 [ - + ][ - - ]: 21863 : Assert( (n[0]*n[0] + n[1]*n[1] + n[2]*n[2] - 1.0) <
[ - - ][ - - ]
463 : : 1.0e+3*std::numeric_limits< tk::real >::epsilon(),
464 : : "Non-unit normal" );
465 : : }
466 : : }
467 : :
468 : : // Replace global->local ids associated to boundary point normals
469 : 459 : decltype(m_bnorm) loc;
470 [ + + ]: 1612 : for (auto& [s,b] : m_bnorm) {
471 [ + - ]: 1153 : auto& bnd = loc[s];
472 [ + + ]: 23016 : for (auto&& [g,n] : b) {
473 [ + - ][ + - ]: 21863 : bnd[ tk::cref_find(lid,g) ] = std::move(n);
474 : : }
475 : : }
476 : 459 : m_bnorm = std::move(loc);
477 : 459 : }
478 : :
479 : : void
480 : 459 : KozCG::streamable()
481 : : // *****************************************************************************
482 : : // Convert integrals into streamable data structures
483 : : // *****************************************************************************
484 : : {
485 : : // Query surface integral output nodes
486 : 459 : std::unordered_map< int, std::vector< std::size_t > > surfintnodes;
487 : 459 : const auto& is = g_cfg.get< tag::integout >();
488 [ + - ]: 459 : std::set< int > outsets( begin(is), end(is) );
489 [ - + ]: 459 : for (auto s : outsets) {
490 [ - - ]: 0 : auto m = m_bface.find(s);
491 [ - - ]: 0 : if (m != end(m_bface)) {
492 [ - - ]: 0 : auto& n = surfintnodes[ m->first ]; // associate set id
493 [ - - ]: 0 : for (auto f : m->second) { // face ids on side set
494 [ - - ]: 0 : n.push_back( m_triinpoel[f*3+0] ); // nodes on side set
495 [ - - ]: 0 : n.push_back( m_triinpoel[f*3+1] );
496 [ - - ]: 0 : n.push_back( m_triinpoel[f*3+2] );
497 : : }
498 : : }
499 : : }
500 [ - - ][ - + ]: 459 : for (auto& [s,n] : surfintnodes) tk::unique( n );
501 : : // Prepare surface integral data
502 : 459 : tk::destroy( m_surfint );
503 [ + - ]: 459 : const auto& gid = Disc()->Gid();
504 [ - + ]: 459 : for (auto&& [s,n] : surfintnodes) {
505 [ - - ]: 0 : auto& sint = m_surfint[s]; // associate set id
506 : 0 : auto& nodes = sint.first;
507 : 0 : auto& ndA = sint.second;
508 : 0 : nodes = std::move(n);
509 [ - - ]: 0 : ndA.resize( nodes.size()*3 );
510 : 0 : std::size_t a = 0;
511 [ - - ]: 0 : for (auto p : nodes) {
512 [ - - ]: 0 : const auto& b = tk::cref_find( m_bndpoinint, gid[p] );
513 : 0 : ndA[a*3+0] = b[0]; // store ni * dA
514 : 0 : ndA[a*3+1] = b[1];
515 : 0 : ndA[a*3+2] = b[2];
516 : 0 : ++a;
517 : : }
518 : : }
519 : :
520 : : // Convert symmetry BC data to streamable data structures
521 : 459 : tk::destroy( m_symbcnodes );
522 : 459 : tk::destroy( m_symbcnorms );
523 [ + + ]: 4699 : for (auto p : m_symbcnodeset) {
524 [ + + ]: 17546 : for (const auto& s : g_cfg.get< tag::bc_sym >()) {
525 [ + - ]: 13306 : auto m = m_bnorm.find(s);
526 [ + - ]: 13306 : if (m != end(m_bnorm)) {
527 [ + - ]: 13306 : auto r = m->second.find(p);
528 [ + + ]: 13306 : if (r != end(m->second)) {
529 [ + - ]: 5268 : m_symbcnodes.push_back( p );
530 [ + - ]: 5268 : m_symbcnorms.push_back( r->second[0] );
531 [ + - ]: 5268 : m_symbcnorms.push_back( r->second[1] );
532 [ + - ]: 5268 : m_symbcnorms.push_back( r->second[2] );
533 : : }
534 : : }
535 : : }
536 : : }
537 : 459 : tk::destroy( m_symbcnodeset );
538 : :
539 : : // Convert farfield BC data to streamable data structures
540 : 459 : tk::destroy( m_farbcnodes );
541 : 459 : tk::destroy( m_farbcnorms );
542 [ - + ]: 459 : for (auto p : m_farbcnodeset) {
543 [ - - ]: 0 : for (const auto& s : g_cfg.get< tag::bc_far >()) {
544 [ - - ]: 0 : auto n = m_bnorm.find(s);
545 [ - - ]: 0 : if (n != end(m_bnorm)) {
546 [ - - ]: 0 : auto a = n->second.find(p);
547 [ - - ]: 0 : if (a != end(n->second)) {
548 [ - - ]: 0 : m_farbcnodes.push_back( p );
549 [ - - ]: 0 : m_farbcnorms.push_back( a->second[0] );
550 [ - - ]: 0 : m_farbcnorms.push_back( a->second[1] );
551 [ - - ]: 0 : m_farbcnorms.push_back( a->second[2] );
552 : : }
553 : : }
554 : : }
555 : : }
556 : 459 : tk::destroy( m_farbcnodeset );
557 : 459 : tk::destroy( m_bnorm );
558 : 459 : }
559 : :
560 : :
561 : : void
562 : : // cppcheck-suppress unusedFunction
563 : 459 : KozCG::merge()
564 : : // *****************************************************************************
565 : : // Combine own and communicated portions of the integrals
566 : : // *****************************************************************************
567 : : {
568 : 459 : auto d = Disc();
569 : :
570 : : // Combine own and communicated contributions to boundary point normals
571 : 459 : bnorm();
572 : :
573 : : // Convert integrals into streamable data structures
574 : 459 : streamable();
575 : :
576 : : // Enforce boundary conditions using (re-)computed boundary data
577 : 459 : BC( m_u, d->T() );
578 : :
579 [ + - ]: 459 : if (d->Initial()) {
580 : : // Output initial conditions to file
581 [ + - ][ + - ]: 459 : writeFields( CkCallback(CkIndex_KozCG::start(), thisProxy[thisIndex]) );
[ + - ][ + - ]
582 : : } else {
583 : 0 : feop_complete();
584 : : }
585 : 459 : }
586 : :
587 : : void
588 : 9465 : KozCG::BC( tk::Fields& u, tk::real t )
589 : : // *****************************************************************************
590 : : // Apply boundary conditions
591 : : //! \param[in,out] u Solution to apply BCs to
592 : : //! \param[in] t Physical time
593 : : // *****************************************************************************
594 : : {
595 : 9465 : auto d = Disc();
596 : :
597 : : // Apply Dirichlet BCs
598 [ + - ]: 9465 : physics::dirbc( u, t, d->Coord(), d->BoxNodes(), m_dirbcmasks );
599 : :
600 : : // Apply symmetry BCs
601 : 9465 : physics::symbc( u, m_symbcnodes, m_symbcnorms, /*pos=*/1 );
602 : :
603 : : // Apply farfield BCs
604 : 9465 : physics::farbc( u, m_farbcnodes, m_farbcnorms );
605 : :
606 : : // Apply pressure BCs
607 : 9465 : physics::prebc( u, m_prebcnodes, m_prebcvals );
608 : 9465 : }
609 : :
610 : : void
611 : 9006 : KozCG::dt()
612 : : // *****************************************************************************
613 : : // Compute time step size
614 : : // *****************************************************************************
615 : : {
616 : 9006 : tk::real mindt = std::numeric_limits< tk::real >::max();
617 : :
618 : 9006 : auto const_dt = g_cfg.get< tag::dt >();
619 : 9006 : auto eps = std::numeric_limits< tk::real >::epsilon();
620 [ + - ]: 9006 : auto d = Disc();
621 : :
622 : : // use constant dt if configured
623 [ + + ]: 9006 : if (std::abs(const_dt) > eps) {
624 : :
625 : : // cppcheck-suppress redundantInitialization
626 : 2920 : mindt = const_dt;
627 : :
628 : : } else {
629 : :
630 : 6086 : const auto& vol = d->Vol();
631 : 6086 : auto cfl = g_cfg.get< tag::cfl >();
632 : :
633 [ - + ]: 6086 : if (g_cfg.get< tag::steady >()) {
634 : :
635 [ - - ]: 0 : for (std::size_t p=0; p<m_u.nunk(); ++p) {
636 [ - - ]: 0 : auto r = m_u(p,0);
637 [ - - ]: 0 : auto u = m_u(p,1)/r;
638 [ - - ]: 0 : auto v = m_u(p,2)/r;
639 [ - - ]: 0 : auto w = m_u(p,3)/r;
640 [ - - ]: 0 : auto pr = eos::pressure( m_u(p,4) - 0.5*r*(u*u + v*v + w*w) );
641 : 0 : auto c = eos::soundspeed( r, std::max(pr,0.0) );
642 : 0 : auto L = std::cbrt( vol[p] );
643 : 0 : auto vel = std::sqrt( u*u + v*v + w*w );
644 : 0 : m_dtp[p] = L / std::max( vel+c, 1.0e-8 ) * cfl;
645 : : }
646 [ - - ]: 0 : mindt = *std::min_element( begin(m_dtp), end(m_dtp) );
647 : :
648 : : } else {
649 : :
650 [ + + ]: 351185 : for (std::size_t p=0; p<m_u.nunk(); ++p) {
651 [ + - ]: 345099 : auto r = m_u(p,0);
652 [ + - ]: 345099 : auto u = m_u(p,1)/r;
653 [ + - ]: 345099 : auto v = m_u(p,2)/r;
654 [ + - ]: 345099 : auto w = m_u(p,3)/r;
655 [ + - ]: 345099 : auto pr = eos::pressure( m_u(p,4) - 0.5*r*(u*u + v*v + w*w) );
656 : 345099 : auto c = eos::soundspeed( r, std::max(pr,0.0) );
657 : 345099 : auto L = std::cbrt( vol[p] );
658 : 345099 : auto vel = std::sqrt( u*u + v*v + w*w );
659 : 345099 : auto euler_dt = L / std::max( vel+c, 1.0e-8 );
660 : 345099 : mindt = std::min( mindt, euler_dt );
661 : : }
662 : 6086 : mindt *= cfl;
663 : :
664 : : }
665 : :
666 : 6086 : mindt *= m_freezeflow;
667 : :
668 : : }
669 : :
670 : : // Actiavate SDAG waits for next time step
671 [ + - ][ + - ]: 9006 : thisProxy[ thisIndex ].wait4rhs();
672 [ + - ][ + - ]: 9006 : thisProxy[ thisIndex ].wait4aec();
673 [ + - ][ + - ]: 9006 : thisProxy[ thisIndex ].wait4alw();
674 [ + - ][ + - ]: 9006 : thisProxy[ thisIndex ].wait4sol();
675 [ + - ][ + - ]: 9006 : thisProxy[ thisIndex ].wait4step();
676 : :
677 : : // Contribute to minimum dt across all chares and advance to next step
678 [ + - ]: 9006 : contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
679 [ + - ][ + - ]: 18012 : CkCallback(CkReductionTarget(KozCG,advance), thisProxy) );
680 : 9006 : }
681 : :
682 : : void
683 : 9006 : KozCG::advance( tk::real newdt )
684 : : // *****************************************************************************
685 : : // Advance equations to next time step
686 : : //! \param[in] newdt The smallest dt across the whole problem
687 : : // *****************************************************************************
688 : : {
689 : : // Set new time step size
690 : 9006 : Disc()->setdt( newdt );
691 : :
692 : : // Compute rhs
693 : 9006 : rhs();
694 : 9006 : }
695 : :
696 : : void
697 : 9006 : KozCG::rhs()
698 : : // *****************************************************************************
699 : : // Compute right-hand side
700 : : // *****************************************************************************
701 : : {
702 : 9006 : auto d = Disc();
703 : :
704 : : // Compute own portion of right-hand side for all equations
705 : 9006 : kozak::rhs( d->Inpoel(), d->Coord(), d->T(), d->Dt(), m_tp, m_dtp, m_u,
706 : 9006 : m_rhs );
707 : :
708 : : // Communicate rhs to other chares on chare-boundary
709 [ + + ]: 9006 : if (d->NodeCommMap().empty()) {
710 : 310 : comrhs_complete();
711 : : } else {
712 : 8696 : const auto& lid = d->Lid();
713 [ + + ]: 91842 : for (const auto& [c,n] : d->NodeCommMap()) {
714 : 83146 : decltype(m_rhsc) exp;
715 [ + - ][ + - ]: 475864 : for (auto g : n) exp[g] = m_rhs[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
716 [ + - ][ + - ]: 83146 : thisProxy[c].comrhs( exp );
717 : 83146 : }
718 : : }
719 : 9006 : ownrhs_complete();
720 : 9006 : }
721 : :
722 : : void
723 : 83146 : KozCG::comrhs(
724 : : const std::unordered_map< std::size_t, std::vector< tk::real > >& inrhs )
725 : : // *****************************************************************************
726 : : // Receive contributions to right-hand side vector on chare-boundaries
727 : : //! \param[in] inrhs Partial contributions of RHS to chare-boundary nodes. Key:
728 : : //! global mesh node IDs, value: contributions for all scalar components.
729 : : // *****************************************************************************
730 : : {
731 : : using tk::operator+=;
732 [ + - ][ + - ]: 475864 : for (const auto& [g,r] : inrhs) m_rhsc[g] += r;
[ + + ]
733 : :
734 : : // When we have heard from all chares we communicate with, this chare is done
735 [ + + ]: 83146 : if (++m_nrhs == Disc()->NodeCommMap().size()) {
736 : 8696 : m_nrhs = 0;
737 : 8696 : comrhs_complete();
738 : : }
739 : 83146 : }
740 : :
741 : : void
742 : 9006 : KozCG::fct()
743 : : // *****************************************************************************
744 : : // Continue with flux-corrected transport if enabled
745 : : // *****************************************************************************
746 : : {
747 : 9006 : auto d = Disc();
748 : 9006 : const auto& lid = d->Lid();
749 : :
750 : : // Combine own and communicated contributions to rhs
751 [ + + ]: 207564 : for (const auto& [g,r] : m_rhsc) {
752 [ + - ]: 198558 : auto i = tk::cref_find( lid, g );
753 [ + - ][ + + ]: 1226488 : for (std::size_t c=0; c<r.size(); ++c) m_rhs(i,c) += r[c];
754 : : }
755 : 9006 : tk::destroy(m_rhsc);
756 : :
757 [ + + ]: 9006 : if (g_cfg.get< tag::fct >()) aec(); else solve();
758 : 9006 : }
759 : :
760 : : void
761 : : // cppcheck-suppress unusedFunction
762 : 4210 : KozCG::aec()
763 : : // *****************************************************************************
764 : : // Compute antidiffusive contributions: P+/-
765 : : // *****************************************************************************
766 : : {
767 : 4210 : auto d = Disc();
768 : 4210 : const auto ncomp = m_u.nprop();
769 : 4210 : const auto& lid = d->Lid();
770 : :
771 : : // Antidiffusive contributions: P+/-
772 : :
773 : 4210 : auto ctau = g_cfg.get< tag::fctdif >();
774 : 4210 : m_p.fill( 0.0 );
775 : :
776 : 4210 : const auto& inpoel = d->Inpoel();
777 : 4210 : const auto& coord = d->Coord();
778 : 4210 : const auto& x = coord[0];
779 : 4210 : const auto& y = coord[1];
780 : 4210 : const auto& z = coord[2];
781 : :
782 [ + + ]: 474910 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
783 : 470700 : const auto N = inpoel.data() + e*4;
784 : : const std::array< tk::real, 3 >
785 : 470700 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
786 : 470700 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
787 : 470700 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
788 : 470700 : const auto J = tk::triple( ba, ca, da );
789 [ + + ]: 3042780 : for (std::size_t c=0; c<ncomp; ++c) {
790 : 2572080 : auto p = c*2;
791 : 2572080 : auto n = p+1;
792 : 2572080 : tk::real aec[4] = { 0.0, 0.0, 0.0, 0.0 };
793 [ + + ]: 12860400 : for (std::size_t a=0; a<4; ++a) {
794 [ + + ]: 51441600 : for (std::size_t b=0; b<4; ++b) {
795 [ + + ]: 41153280 : auto m = J/120.0 * ((a == b) ? 3.0 : -1.0);
796 [ + - ]: 41153280 : aec[a] += m * ctau * m_u(N[b],c);
797 : : }
798 [ + - ]: 10288320 : m_p(N[a],p) += std::max(0.0,aec[a]);
799 [ + - ]: 10288320 : m_p(N[a],n) += std::min(0.0,aec[a]);
800 : : }
801 : : }
802 : : }
803 : :
804 : : // Apply symmetry BCs on AEC
805 [ + + ]: 56890 : for (std::size_t i=0; i<m_symbcnodes.size(); ++i) {
806 : 52680 : auto p = m_symbcnodes[i];
807 : 52680 : auto nx = m_symbcnorms[i*3+0];
808 : 52680 : auto ny = m_symbcnorms[i*3+1];
809 : 52680 : auto nz = m_symbcnorms[i*3+2];
810 : 52680 : auto rvnp = m_p(p,2)*nx + m_p(p,4)*ny + m_p(p,6)*nz;
811 : 52680 : auto rvnn = m_p(p,3)*nx + m_p(p,5)*ny + m_p(p,7)*nz;
812 : 52680 : m_p(p,2) -= rvnp * nx;
813 : 52680 : m_p(p,3) -= rvnn * nx;
814 : 52680 : m_p(p,4) -= rvnp * ny;
815 : 52680 : m_p(p,5) -= rvnn * ny;
816 : 52680 : m_p(p,6) -= rvnp * nz;
817 : 52680 : m_p(p,7) -= rvnn * nz;
818 : : }
819 : :
820 : : // Communicate antidiffusive edge and low-order solution contributions
821 [ + + ]: 4210 : if (d->NodeCommMap().empty()) {
822 : 100 : comaec_complete();
823 : : } else {
824 [ + + ]: 44410 : for (const auto& [c,n] : d->NodeCommMap()) {
825 : 40300 : decltype(m_pc) exp;
826 [ + - ][ + - ]: 223540 : for (auto g : n) exp[g] = m_p[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
827 [ + - ][ + - ]: 40300 : thisProxy[c].comaec( exp );
828 : 40300 : }
829 : : }
830 : 4210 : ownaec_complete();
831 : 4210 : }
832 : :
833 : : void
834 : 40300 : KozCG::comaec( const std::unordered_map< std::size_t,
835 : : std::vector< tk::real > >& inaec )
836 : : // *****************************************************************************
837 : : // Receive antidiffusive and low-order contributions on chare-boundaries
838 : : //! \param[in] inaec Partial contributions of antidiffusive edge and low-order
839 : : //! solution contributions on chare-boundary nodes. Key: global mesh node IDs,
840 : : //! value: 0: antidiffusive contributions, 1: low-order solution.
841 : : // *****************************************************************************
842 : : {
843 : : using tk::operator+=;
844 [ + - ][ + - ]: 223540 : for (const auto& [g,a] : inaec) m_pc[g] += a;
[ + + ]
845 : :
846 : : // When we have heard from all chares we communicate with, this chare is done
847 [ + + ]: 40300 : if (++m_naec == Disc()->NodeCommMap().size()) {
848 : 4110 : m_naec = 0;
849 : 4110 : comaec_complete();
850 : : }
851 : 40300 : }
852 : :
853 : : void
854 : 4210 : KozCG::alw()
855 : : // *****************************************************************************
856 : : // Compute allowed limits, Q+/-
857 : : // *****************************************************************************
858 : : {
859 : 4210 : auto d = Disc();
860 : 4210 : const auto steady = g_cfg.get< tag::steady >();
861 : 4210 : const auto npoin = m_u.nunk();
862 : 4210 : const auto ncomp = m_u.nprop();
863 : 4210 : const auto& lid = d->Lid();
864 : 4210 : const auto& vol = d->Vol();
865 : 4210 : const auto& inpoel = d->Inpoel();
866 : :
867 : : // Combine own and communicated contributions to antidiffusive contributions
868 : : // and low-order solution
869 [ + + ]: 99880 : for (const auto& [g,p] : m_pc) {
870 [ + - ]: 95670 : auto i = tk::cref_find( lid, g );
871 [ + - ][ + + ]: 1122650 : for (std::size_t c=0; c<p.size(); ++c) m_p(i,c) += p[c];
872 : : }
873 : 4210 : tk::destroy(m_pc);
874 : :
875 : : // Finish computing antidiffusive contributions and low-order solution
876 : 4210 : auto dt = d->Dt();
877 [ + + ]: 215140 : for (std::size_t i=0; i<npoin; ++i) {
878 [ - + ]: 210930 : if (steady) dt = m_dtp[i];
879 [ + + ]: 1356600 : for (std::size_t c=0; c<ncomp; ++c) {
880 : 1145670 : auto p = c*2;
881 : 1145670 : auto n = p+1;
882 : 1145670 : m_p(i,p) /= vol[i];
883 : 1145670 : m_p(i,n) /= vol[i];
884 : : // low-order solution
885 : 1145670 : m_rhs(i,c) = m_u(i,c) + dt*m_rhs(i,c)/vol[i] - m_p(i,p) - m_p(i,n);
886 : : }
887 : : }
888 : :
889 : : // Allowed limits: Q+/-
890 : :
891 : : using std::max;
892 : : using std::min;
893 : :
894 : 4210 : auto large = std::numeric_limits< tk::real >::max();
895 [ + + ]: 215140 : for (std::size_t i=0; i<m_q.nunk(); ++i) {
896 [ + + ]: 1356600 : for (std::size_t c=0; c<m_q.nprop()/2; ++c) {
897 : 1145670 : m_q(i,c*2+0) = -large;
898 : 1145670 : m_q(i,c*2+1) = +large;
899 : : }
900 : : }
901 : :
902 [ + + ]: 474910 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
903 : 470700 : const auto N = inpoel.data() + e*4;
904 [ + + ]: 3042780 : for (std::size_t c=0; c<ncomp; ++c) {
905 : 2572080 : auto alwp = -large;
906 : 2572080 : auto alwn = +large;
907 [ + + ]: 12860400 : for (std::size_t a=0; a<4; ++a) {
908 [ + + ]: 10288320 : if (g_cfg.get< tag::fctclip >()) {
909 [ + - ]: 2122400 : alwp = max( alwp, m_rhs(N[a],c) );
910 [ + - ]: 2122400 : alwn = min( alwn, m_rhs(N[a],c) );
911 : : } else {
912 [ + - ][ + - ]: 8165920 : alwp = max( alwp, max(m_rhs(N[a],c), m_u(N[a],c)) );
913 [ + - ][ + - ]: 8165920 : alwn = min( alwn, min(m_rhs(N[a],c), m_u(N[a],c)) );
914 : : }
915 : : }
916 : 2572080 : auto p = c*2;
917 : 2572080 : auto n = p+1;
918 [ + + ]: 12860400 : for (std::size_t a=0; a<4; ++a) {
919 [ + - ][ + - ]: 10288320 : m_q(N[a],p) = max(m_q(N[a],p), alwp);
920 [ + - ][ + - ]: 10288320 : m_q(N[a],n) = min(m_q(N[a],n), alwn);
921 : : }
922 : : }
923 : : }
924 : :
925 : : // Communicate allowed limits contributions
926 [ + + ]: 4210 : if (d->NodeCommMap().empty()) {
927 : 100 : comalw_complete();
928 : : } else {
929 [ + + ]: 44410 : for (const auto& [c,n] : d->NodeCommMap()) {
930 : 40300 : decltype(m_qc) exp;
931 [ + - ][ + - ]: 223540 : for (auto g : n) exp[g] = m_q[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
932 [ + - ][ + - ]: 40300 : thisProxy[c].comalw( exp );
933 : 40300 : }
934 : : }
935 : 4210 : ownalw_complete();
936 : 4210 : }
937 : :
938 : : void
939 : 40300 : KozCG::comalw( const std::unordered_map< std::size_t,
940 : : std::vector< tk::real > >& inalw )
941 : : // *****************************************************************************
942 : : // Receive allowed limits contributions on chare-boundaries
943 : : //! \param[in] inalw Partial contributions of allowed limits contributions on
944 : : //! chare-boundary nodes. Key: global mesh node IDs, value: allowed limit
945 : : //! contributions.
946 : : // *****************************************************************************
947 : : {
948 [ + + ]: 223540 : for (const auto& [g,alw] : inalw) {
949 [ + - ]: 183240 : auto& q = m_qc[g];
950 [ + - ]: 183240 : q.resize( alw.size() );
951 [ + + ]: 1143840 : for (std::size_t c=0; c<alw.size()/2; ++c) {
952 : 960600 : auto p = c*2;
953 : 960600 : auto n = p+1;
954 : 960600 : q[p] = std::max( q[p], alw[p] );
955 : 960600 : q[n] = std::min( q[n], alw[n] );
956 : : }
957 : : }
958 : :
959 : : // When we have heard from all chares we communicate with, this chare is done
960 [ + + ]: 40300 : if (++m_nalw == Disc()->NodeCommMap().size()) {
961 : 4110 : m_nalw = 0;
962 : 4110 : comalw_complete();
963 : : }
964 : 40300 : }
965 : :
966 : : void
967 : 4210 : KozCG::lim()
968 : : // *****************************************************************************
969 : : // Compute limit coefficients
970 : : // *****************************************************************************
971 : : {
972 [ + - ]: 4210 : auto d = Disc();
973 : 4210 : const auto npoin = m_u.nunk();
974 : 4210 : const auto ncomp = m_u.nprop();
975 : 4210 : const auto& lid = d->Lid();
976 : :
977 : : using std::max;
978 : : using std::min;
979 : :
980 : : // Combine own and communicated contributions to allowed limits
981 [ + + ]: 99880 : for (const auto& [g,alw] : m_qc) {
982 [ + - ]: 95670 : auto i = tk::cref_find( lid, g );
983 [ + + ]: 609160 : for (std::size_t c=0; c<alw.size()/2; ++c) {
984 : 513490 : auto p = c*2;
985 : 513490 : auto n = p+1;
986 [ + - ][ + - ]: 513490 : m_q(i,p) = max( m_q(i,p), alw[p] );
987 [ + - ][ + - ]: 513490 : m_q(i,n) = min( m_q(i,n), alw[n] );
988 : : }
989 : : }
990 : 4210 : tk::destroy(m_qc);
991 : :
992 : : // Finish computing allowed limits
993 [ + + ]: 215140 : for (std::size_t i=0; i<npoin; ++i) {
994 [ + + ]: 1356600 : for (std::size_t c=0; c<ncomp; ++c) {
995 : 1145670 : auto p = c*2;
996 : 1145670 : auto n = p+1;
997 [ + - ][ + - ]: 1145670 : m_q(i,p) -= m_rhs(i,c);
998 [ + - ][ + - ]: 1145670 : m_q(i,n) -= m_rhs(i,c);
999 : : }
1000 : : }
1001 : :
1002 : : // Limit coefficients, C
1003 : :
1004 [ + + ]: 215140 : for (std::size_t i=0; i<npoin; ++i) {
1005 [ + + ]: 1356600 : for (std::size_t c=0; c<ncomp; ++c) {
1006 : 1145670 : auto p = c*2;
1007 : 1145670 : auto n = p+1;
1008 : 1145670 : auto eps = std::numeric_limits< tk::real >::epsilon();
1009 [ + - ][ + + ]: 1145670 : m_q(i,p) = m_p(i,p) < eps ? 0.0 : min(1.0, m_q(i,p)/m_p(i,p));
[ + - ][ + - ]
[ + - ]
1010 [ + - ][ + + ]: 1145670 : m_q(i,n) = m_p(i,n) > -eps ? 0.0 : min(1.0, m_q(i,n)/m_p(i,n));
[ + - ][ + - ]
[ + - ]
1011 : : }
1012 : : }
1013 : :
1014 : : // Limited antidiffusive contributions
1015 : :
1016 : 4210 : auto ctau = g_cfg.get< tag::fctdif >();
1017 [ + - ]: 4210 : m_a.fill( 0.0 );
1018 : :
1019 : 4210 : const auto& inpoel = d->Inpoel();
1020 : 4210 : const auto& coord = d->Coord();
1021 : 4210 : const auto& x = coord[0];
1022 : 4210 : const auto& y = coord[1];
1023 : 4210 : const auto& z = coord[2];
1024 : :
1025 [ + - ]: 4210 : auto fctsys = g_cfg.get< tag::fctsys >();
1026 [ + + ]: 4990 : for (auto& c : fctsys) --c;
1027 : :
1028 : : #if defined(__clang__)
1029 : : #pragma clang diagnostic push
1030 : : #pragma clang diagnostic ignored "-Wvla"
1031 : : #pragma clang diagnostic ignored "-Wvla-extension"
1032 : : #elif defined(STRICT_GNUC)
1033 : : #pragma GCC diagnostic push
1034 : : #pragma GCC diagnostic ignored "-Wvla"
1035 : : #endif
1036 : :
1037 [ + + ]: 474910 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
1038 : 941400 : const auto N = inpoel.data() + e*4;
1039 : : const std::array< tk::real, 3 >
1040 : 470700 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
1041 : 470700 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
1042 : 470700 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
1043 : 470700 : const auto J = tk::triple( ba, ca, da );
1044 : 470700 : tk::real coef[ncomp], aec[ncomp][4];
1045 [ + + ]: 3042780 : for (std::size_t c=0; c<ncomp; ++c) {
1046 : 2572080 : auto p = c*2;
1047 : 2572080 : auto n = p+1;
1048 : 2572080 : coef[c] = 1.0;
1049 [ + + ]: 12860400 : for (std::size_t a=0; a<4; ++a) {
1050 : 10288320 : aec[c][a] = 0.0;
1051 [ + + ]: 51441600 : for (std::size_t b=0; b<4; ++b) {
1052 [ + + ]: 41153280 : auto m = J/120.0 * ((a == b) ? 3.0 : -1.0);
1053 [ + - ]: 41153280 : aec[c][a] += m * ctau * m_u(N[b],c);
1054 : : }
1055 [ + + ][ + - ]: 10288320 : coef[c] = min(coef[c], aec[c][a] > 0.0 ? m_q(N[a],p) : m_q(N[a],n));
[ + - ]
1056 : : }
1057 : : }
1058 : 470700 : tk::real cs = 1.0;
1059 [ + + ]: 789060 : for (auto c : fctsys) cs = min( cs, coef[c] );
1060 [ + + ]: 789060 : for (auto c : fctsys) coef[c] = cs;
1061 [ + + ]: 3042780 : for (std::size_t c=0; c<ncomp; ++c) {
1062 [ + + ]: 12860400 : for (std::size_t a=0; a<4; ++a) {
1063 [ + - ]: 10288320 : m_a(N[a],c) += coef[c] * aec[c][a];
1064 : : }
1065 : : }
1066 : 470700 : }
1067 : :
1068 : : #if defined(__clang__)
1069 : : #pragma clang diagnostic pop
1070 : : #elif defined(STRICT_GNUC)
1071 : : #pragma GCC diagnostic pop
1072 : : #endif
1073 : :
1074 : : // Communicate limited antidiffusive contributions
1075 [ + + ]: 4210 : if (d->NodeCommMap().empty()) {
1076 [ + - ]: 100 : comlim_complete();
1077 : : } else {
1078 [ + + ]: 44410 : for (const auto& [c,n] : d->NodeCommMap()) {
1079 : 40300 : decltype(m_ac) exp;
1080 [ + - ][ + - ]: 223540 : for (auto g : n) exp[g] = m_a[ tk::cref_find(lid,g) ];
[ + - ][ + + ]
1081 [ + - ][ + - ]: 40300 : thisProxy[c].comlim( exp );
1082 : 40300 : }
1083 : : }
1084 [ + - ]: 4210 : ownlim_complete();
1085 : 4210 : }
1086 : :
1087 : : void
1088 : 40300 : KozCG::comlim( const std::unordered_map< std::size_t,
1089 : : std::vector< tk::real > >& inlim )
1090 : : // *****************************************************************************
1091 : : // Receive limited antidiffusive contributions on chare-boundaries
1092 : : //! \param[in] inlim Partial contributions of limited contributions on
1093 : : //! chare-boundary nodes. Key: global mesh node IDs, value: limited
1094 : : //! contributions.
1095 : : // *****************************************************************************
1096 : : {
1097 : : using tk::operator+=;
1098 [ + - ][ + - ]: 223540 : for (const auto& [g,a] : inlim) m_ac[g] += a;
[ + + ]
1099 : :
1100 : : // When we have heard from all chares we communicate with, this chare is done
1101 [ + + ]: 40300 : if (++m_nlim == Disc()->NodeCommMap().size()) {
1102 : 4110 : m_nlim = 0;
1103 : 4110 : comlim_complete();
1104 : : }
1105 : 40300 : }
1106 : :
1107 : : void
1108 : 9006 : KozCG::solve()
1109 : : // *****************************************************************************
1110 : : // Compute limit coefficients
1111 : : // *****************************************************************************
1112 : : {
1113 [ + - ]: 9006 : auto d = Disc();
1114 : 9006 : const auto npoin = m_u.nunk();
1115 : 9006 : const auto ncomp = m_u.nprop();
1116 : 9006 : const auto& lid = d->Lid();
1117 : 9006 : const auto& vol = d->Vol();
1118 : 9006 : const auto steady = g_cfg.get< tag::steady >();
1119 : :
1120 : : // Combine own and communicated contributions to limited antidiffusive
1121 : : // contributions
1122 [ + + ]: 104676 : for (const auto& [g,a] : m_ac) {
1123 [ + - ]: 95670 : auto i = tk::cref_find( lid, g );
1124 [ + - ][ + + ]: 609160 : for (std::size_t c=0; c<a.size(); ++c) m_a(i,c) += a[c];
1125 : : }
1126 : 9006 : tk::destroy(m_ac);
1127 : :
1128 : 9006 : tk::Fields u;
1129 [ + + ]: 9006 : std::size_t cstart = m_freezeflow > 1.0 ? 5 : 0;
1130 [ + + ][ + - ]: 9006 : if (cstart) u = m_u;
1131 : :
1132 [ + + ]: 9006 : if (g_cfg.get< tag::fct >()) {
1133 : : // Apply limited antidiffusive contributions to low-order solution
1134 [ + + ]: 215140 : for (std::size_t i=0; i<npoin; ++i) {
1135 [ + + ]: 1356600 : for (std::size_t c=0; c<ncomp; ++c) {
1136 [ + - ][ + - ]: 1145670 : m_a(i,c) = m_rhs(i,c) + m_a(i,c)/vol[i];
[ + - ]
1137 : : }
1138 : : }
1139 : : } else {
1140 : : // Apply rhs
1141 : 4796 : auto dt = d->Dt();
1142 [ + + ]: 215945 : for (std::size_t i=0; i<npoin; ++i) {
1143 [ - + ]: 211149 : if (steady) dt = m_dtp[i];
1144 [ + + ]: 1266894 : for (std::size_t c=0; c<ncomp; ++c) {
1145 [ + - ][ + - ]: 1055745 : m_a(i,c) = m_u(i,c) + dt*m_rhs(i,c)/vol[i];
[ + - ]
1146 : : }
1147 : : }
1148 : : }
1149 : :
1150 : : // Apply scalar source to solution (if defined)
1151 [ + - ]: 9006 : auto src = problems::PHYS_SRC();
1152 [ - + ][ - - ]: 9006 : if (src) src( d->Coord(), d->T(), m_a );
1153 : :
1154 : : // Freeze flow if configured and apply multiplier on scalar(s)
1155 [ + + ]: 9006 : if (d->T() > g_cfg.get< tag::freezetime >()) {
1156 : 8547 : m_freezeflow = g_cfg.get< tag::freezeflow >();
1157 : : }
1158 : :
1159 : : // Enforce boundary conditions
1160 [ + - ]: 9006 : BC( m_a, d->T() + d->Dt() );
1161 : :
1162 : : // Explicitly zero out flow for freezeflow
1163 [ + + ]: 9006 : if (cstart) {
1164 [ + + ]: 82800 : for (std::size_t i=0; i<npoin; ++i) {
1165 [ + + ]: 491508 : for (std::size_t c=0; c<cstart; ++c) {
1166 [ + - ][ + - ]: 409590 : m_a(i,c) = u(i,c);
1167 : : }
1168 : : }
1169 : : }
1170 : :
1171 : : // Compute diagnostics, e.g., residuals
1172 : 9006 : auto diag_iter = g_cfg.get< tag::diag_iter >();
1173 [ + - ]: 9006 : auto diag = m_diag.rhocompute( *d, m_a, m_u, diag_iter );
1174 : :
1175 : : // Update solution
1176 [ + - ]: 9006 : m_u = m_a;
1177 [ + - ]: 9006 : m_a.fill( 0.0 );
1178 : :
1179 : : // Increase number of iterations and physical time
1180 [ + - ]: 9006 : d->next();
1181 : :
1182 : : // Advance physical time for local time stepping
1183 [ - + ]: 9006 : if (steady) {
1184 : : using tk::operator+=;
1185 [ - - ]: 0 : m_tp += m_dtp;
1186 : : }
1187 : :
1188 : : // Evaluate residuals
1189 [ + + ][ + - ]: 9006 : if (!diag) evalres( std::vector< tk::real >( ncomp, 1.0 ) );
[ + - ]
1190 : 9006 : }
1191 : :
1192 : : void
1193 : 9006 : KozCG::evalres( const std::vector< tk::real >& l2res )
1194 : : // *****************************************************************************
1195 : : // Evaluate residuals
1196 : : //! \param[in] l2res L2-norms of the residual for each scalar component
1197 : : //! computed across the whole problem
1198 : : // *****************************************************************************
1199 : : {
1200 [ - + ]: 9006 : if (g_cfg.get< tag::steady >()) {
1201 : 0 : const auto rc = g_cfg.get< tag::rescomp >() - 1;
1202 : 0 : Disc()->residual( l2res[rc] );
1203 : : }
1204 : :
1205 : 9006 : refine();
1206 : 9006 : }
1207 : :
1208 : : void
1209 : 9006 : KozCG::refine()
1210 : : // *****************************************************************************
1211 : : // Optionally refine/derefine mesh
1212 : : // *****************************************************************************
1213 : : {
1214 : 9006 : auto d = Disc();
1215 : :
1216 : : // See if this is the last time step
1217 [ + + ]: 9006 : if (d->finished()) m_finished = 1;
1218 : :
1219 : 9006 : auto dtref = g_cfg.get< tag::href_dt >();
1220 : 9006 : auto dtfreq = g_cfg.get< tag::href_dtfreq >();
1221 : :
1222 : : // if t>0 refinement enabled and we hit the frequency
1223 [ - + ][ - - ]: 9006 : if (dtref && !(d->It() % dtfreq)) { // refine
[ - + ]
1224 : :
1225 : 0 : d->refined() = 1;
1226 : 0 : d->startvol();
1227 : 0 : d->Ref()->dtref( m_bface, m_bnode, m_triinpoel );
1228 : :
1229 : : // Activate SDAG waits for re-computing the integrals
1230 [ - - ][ - - ]: 0 : thisProxy[ thisIndex ].wait4int();
1231 : :
1232 : : } else { // do not refine
1233 : :
1234 : 9006 : d->refined() = 0;
1235 : 9006 : feop_complete();
1236 : 9006 : resize_complete();
1237 : :
1238 : : }
1239 : 9006 : }
1240 : :
1241 : : void
1242 : 0 : KozCG::resizePostAMR(
1243 : : const std::vector< std::size_t >& /*ginpoel*/,
1244 : : const tk::UnsMesh::Chunk& chunk,
1245 : : const tk::UnsMesh::Coords& coord,
1246 : : const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& addedNodes,
1247 : : const std::unordered_map< std::size_t, std::size_t >& /*addedTets*/,
1248 : : const std::set< std::size_t >& removedNodes,
1249 : : const std::unordered_map< int, std::unordered_set< std::size_t > >&
1250 : : nodeCommMap,
1251 : : const std::map< int, std::vector< std::size_t > >& bface,
1252 : : const std::map< int, std::vector< std::size_t > >& bnode,
1253 : : const std::vector< std::size_t >& triinpoel )
1254 : : // *****************************************************************************
1255 : : // Receive new mesh from Refiner
1256 : : //! \param[in] ginpoel Mesh connectivity with global node ids
1257 : : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
1258 : : //! \param[in] coord New mesh node coordinates
1259 : : //! \param[in] addedNodes Newly added mesh nodes and their parents (local ids)
1260 : : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
1261 : : //! \param[in] removedNodes Newly removed mesh node local ids
1262 : : //! \param[in] nodeCommMap New node communication map
1263 : : //! \param[in] bface Boundary-faces mapped to side set ids
1264 : : //! \param[in] bnode Boundary-node lists mapped to side set ids
1265 : : //! \param[in] triinpoel Boundary-face connectivity
1266 : : // *****************************************************************************
1267 : : {
1268 [ - - ]: 0 : auto d = Disc();
1269 : :
1270 : 0 : d->Itf() = 0; // Zero field output iteration count if AMR
1271 : 0 : ++d->Itr(); // Increase number of iterations with a change in the mesh
1272 : :
1273 : : // Resize mesh data structures after mesh refinement
1274 [ - - ]: 0 : d->resizePostAMR( chunk, coord, nodeCommMap, removedNodes );
1275 : :
1276 [ - - ][ - - ]: 0 : Assert(coord[0].size() == m_u.nunk()-removedNodes.size()+addedNodes.size(),
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
1277 : : "Incorrect vector length post-AMR: expected length after resizing = " +
1278 : : std::to_string(coord[0].size()) + ", actual unknown vector length = " +
1279 : : std::to_string(m_u.nunk()-removedNodes.size()+addedNodes.size()));
1280 : :
1281 : : // Remove newly removed nodes from solution vectors
1282 [ - - ]: 0 : m_u.rm( removedNodes );
1283 [ - - ]: 0 : m_rhs.rm( removedNodes );
1284 : :
1285 : : // Resize auxiliary solution vectors
1286 : 0 : auto npoin = coord[0].size();
1287 [ - - ]: 0 : m_u.resize( npoin );
1288 [ - - ]: 0 : m_rhs.resize( npoin );
1289 : :
1290 : : // Update solution on new mesh
1291 [ - - ]: 0 : for (const auto& n : addedNodes)
1292 [ - - ]: 0 : for (std::size_t c=0; c<m_u.nprop(); ++c) {
1293 [ - - ][ - - ]: 0 : Assert(n.first < m_u.nunk(), "Added node index out of bounds post-AMR");
[ - - ][ - - ]
1294 [ - - ][ - - ]: 0 : Assert(n.second[0] < m_u.nunk() && n.second[1] < m_u.nunk(),
[ - - ][ - - ]
[ - - ]
1295 : : "Indices of parent-edge nodes out of bounds post-AMR");
1296 [ - - ][ - - ]: 0 : m_u(n.first,c) = (m_u(n.second[0],c) + m_u(n.second[1],c))/2.0;
[ - - ]
1297 : : }
1298 : :
1299 : : // Update physical-boundary node-, face-, and element lists
1300 [ - - ]: 0 : m_bnode = bnode;
1301 [ - - ]: 0 : m_bface = bface;
1302 [ - - ]: 0 : m_triinpoel = tk::remap( triinpoel, d->Lid() );
1303 : :
1304 : 0 : auto meshid = d->MeshId();
1305 [ - - ]: 0 : contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1306 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,resized), d->Tr()) );
1307 : 0 : }
1308 : :
1309 : : void
1310 : 1402 : KozCG::writeFields( CkCallback cb )
1311 : : // *****************************************************************************
1312 : : // Output mesh-based fields to file
1313 : : //! \param[in] cb Function to continue with after the write
1314 : : // *****************************************************************************
1315 : : {
1316 [ + + ][ + - ]: 1402 : if (g_cfg.get< tag::benchmark >()) { cb.send(); return; }
1317 : :
1318 [ + - ]: 818 : auto d = Disc();
1319 : 818 : auto ncomp = m_u.nprop();
1320 : :
1321 : : // Field output
1322 : :
1323 : : std::vector< std::string > nodefieldnames
1324 [ + - ][ + + ]: 5726 : {"density", "velocityx", "velocityy", "velocityz", "energy", "pressure"};
[ - - ]
1325 [ - + ][ - - ]: 818 : if (g_cfg.get< tag::steady >()) nodefieldnames.push_back( "mach" );
[ - - ]
1326 : :
1327 : : using tk::operator/=;
1328 [ + - ]: 818 : auto r = m_u.extract(0);
1329 [ + - ][ + - ]: 818 : auto u = m_u.extract(1); u /= r;
1330 [ + - ][ + - ]: 818 : auto v = m_u.extract(2); v /= r;
1331 [ + - ][ + - ]: 818 : auto w = m_u.extract(3); w /= r;
1332 [ + - ][ + - ]: 818 : auto e = m_u.extract(4); e /= r;
1333 [ + - ]: 818 : std::vector< tk::real > pr( m_u.nunk() ), ma;
1334 [ - + ][ - - ]: 818 : if (g_cfg.get< tag::steady >()) ma.resize( m_u.nunk() );
1335 [ + + ]: 56795 : for (std::size_t i=0; i<pr.size(); ++i) {
1336 : 55977 : auto vv = u[i]*u[i] + v[i]*v[i] + w[i]*w[i];
1337 : 55977 : pr[i] = eos::pressure( r[i]*(e[i] - 0.5*vv) );
1338 [ - + ]: 55977 : if (g_cfg.get< tag::steady >()) {
1339 : 0 : ma[i] = std::sqrt(vv) / eos::soundspeed( r[i], pr[i] );
1340 : : }
1341 : : }
1342 : :
1343 : : std::vector< std::vector< tk::real > > nodefields{
1344 : 4090 : std::move(r), std::move(u), std::move(v), std::move(w), std::move(e),
1345 [ + - ][ + + ]: 5726 : std::move(pr) };
[ - - ]
1346 [ - + ][ - - ]: 818 : if (g_cfg.get< tag::steady >()) nodefields.push_back( std::move(ma) );
1347 : :
1348 [ + + ]: 916 : for (std::size_t c=0; c<ncomp-5; ++c) {
1349 [ + - ][ + - ]: 98 : nodefieldnames.push_back( "c" + std::to_string(c) );
[ + - ]
1350 [ + - ][ + - ]: 98 : nodefields.push_back( m_u.extract(5+c) );
1351 : : }
1352 : :
1353 : : // query function to evaluate analytic solution (if defined)
1354 [ + - ]: 818 : auto sol = problems::SOL();
1355 : :
1356 [ + + ]: 818 : if (sol) {
1357 : 755 : const auto& coord = d->Coord();
1358 : 755 : const auto& x = coord[0];
1359 : 755 : const auto& y = coord[1];
1360 : 755 : const auto& z = coord[2];
1361 [ + - ]: 755 : auto an = m_u;
1362 [ + - ]: 755 : std::vector< tk::real > ap( m_u.nunk() );
1363 [ + + ]: 45763 : for (std::size_t i=0; i<an.nunk(); ++i) {
1364 [ + - ]: 45008 : auto s = sol( x[i], y[i], z[i], d->T() );
1365 : 45008 : s[1] /= s[0];
1366 : 45008 : s[2] /= s[0];
1367 : 45008 : s[3] /= s[0];
1368 : 45008 : s[4] /= s[0];
1369 [ + - ][ + + ]: 279150 : for (std::size_t c=0; c<s.size(); ++c) an(i,c) = s[c];
1370 : 45008 : s[4] -= 0.5*(s[1]*s[1] + s[2]*s[2] + s[3]*s[3]);
1371 : 45008 : ap[i] = eos::pressure( s[0]*s[4] );
1372 : 45008 : }
1373 [ + + ]: 4530 : for (std::size_t c=0; c<5; ++c) {
1374 [ + - ][ + - ]: 3775 : nodefieldnames.push_back( nodefieldnames[c] + "_analytic" );
1375 [ + - ][ + - ]: 3775 : nodefields.push_back( an.extract(c) );
1376 : : }
1377 [ + - ][ + - ]: 755 : nodefieldnames.push_back( nodefieldnames[5] + "_analytic" );
1378 [ + - ]: 755 : nodefields.push_back( std::move(ap) );
1379 [ + + ]: 853 : for (std::size_t c=0; c<ncomp-5; ++c) {
1380 [ + - ][ + - ]: 98 : nodefieldnames.push_back( nodefieldnames[6+c] + "_analytic" );
1381 [ + - ][ + - ]: 98 : nodefields.push_back( an.extract(5+c) );
1382 : : }
1383 : 755 : }
1384 : :
1385 [ - + ][ - - ]: 818 : Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
[ - - ][ - - ]
1386 : :
1387 : : // Surface output
1388 : :
1389 : 818 : std::vector< std::string > nodesurfnames;
1390 : 818 : std::vector< std::vector< tk::real > > nodesurfs;
1391 : :
1392 : 818 : const auto& f = g_cfg.get< tag::fieldout >();
1393 : :
1394 [ + + ]: 818 : if (!f.empty()) {
1395 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "density" );
1396 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "velocityx" );
1397 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "velocityy" );
1398 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "velocityz" );
1399 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "energy" );
1400 [ + - ][ + - ]: 19 : nodesurfnames.push_back( "pressure" );
1401 : :
1402 [ - + ]: 19 : for (std::size_t c=0; c<ncomp-5; ++c) {
1403 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "c" + std::to_string(c) );
[ - - ]
1404 : : }
1405 [ - + ]: 19 : if (g_cfg.get< tag::steady >()) {
1406 [ - - ][ - - ]: 0 : nodesurfnames.push_back( "mach" );
1407 : : }
1408 : :
1409 [ + - ]: 19 : auto bnode = tk::bfacenodes( m_bface, m_triinpoel );
1410 [ + - ]: 19 : std::set< int > outsets( begin(f), end(f) );
1411 [ + + ]: 71 : for (auto sideset : outsets) {
1412 [ + - ]: 52 : auto b = bnode.find(sideset);
1413 [ + + ]: 52 : if (b == end(bnode)) continue;
1414 : 46 : const auto& nodes = b->second;
1415 : 46 : auto i = nodesurfs.size();
1416 : 46 : auto ns = ncomp + 1;
1417 [ - + ]: 46 : if (g_cfg.get< tag::steady >()) ++ns;
1418 [ + - ]: 46 : nodesurfs.insert( end(nodesurfs), ns,
1419 [ + - ]: 92 : std::vector< tk::real >( nodes.size() ) );
1420 : 46 : std::size_t j = 0;
1421 [ + + ]: 4876 : for (auto n : nodes) {
1422 [ + - ]: 4830 : const auto s = m_u[n];
1423 : 4830 : std::size_t p = 0;
1424 : 4830 : nodesurfs[i+(p++)][j] = s[0];
1425 : 4830 : nodesurfs[i+(p++)][j] = s[1]/s[0];
1426 : 4830 : nodesurfs[i+(p++)][j] = s[2]/s[0];
1427 : 4830 : nodesurfs[i+(p++)][j] = s[3]/s[0];
1428 : 4830 : nodesurfs[i+(p++)][j] = s[4]/s[0];
1429 : 4830 : auto vv = (s[1]*s[1] + s[2]*s[2] + s[3]*s[3])/s[0]/s[0];
1430 : 4830 : auto ei = s[4]/s[0] - 0.5*vv;
1431 : 4830 : auto sp = eos::pressure( s[0]*ei );
1432 : 4830 : nodesurfs[i+(p++)][j] = sp;
1433 [ - + ]: 4830 : for (std::size_t c=0; c<ncomp-5; ++c) nodesurfs[i+(p++)+c][j] = s[5+c];
1434 [ - + ]: 4830 : if (g_cfg.get< tag::steady >()) {
1435 : 0 : nodesurfs[i+(p++)][j] = std::sqrt(vv) / eos::soundspeed( s[0], sp );
1436 : : }
1437 : 4830 : ++j;
1438 : 4830 : }
1439 : : }
1440 : 19 : }
1441 : :
1442 : : // Send mesh and fields data (solution dump) for output to file
1443 [ + - ][ + - ]: 1636 : d->write( d->Inpoel(), d->Coord(), m_bface, tk::remap(m_bnode,d->Lid()),
1444 : 818 : m_triinpoel, {}, nodefieldnames, {}, nodesurfnames,
1445 : : {}, nodefields, {}, nodesurfs, cb );
1446 [ + - ][ + - ]: 3272 : }
[ + - ][ + - ]
[ + - ][ + - ]
[ - - ][ - - ]
[ - - ][ - - ]
1447 : :
1448 : : void
1449 : 9006 : KozCG::out()
1450 : : // *****************************************************************************
1451 : : // Output mesh field data
1452 : : // *****************************************************************************
1453 : : {
1454 : 9006 : auto d = Disc();
1455 : :
1456 : : // Time history
1457 [ + + ][ + + ]: 9006 : if (d->histiter() or d->histtime() or d->histrange()) {
[ - + ][ + + ]
1458 : 214 : auto ncomp = m_u.nprop();
1459 : 214 : const auto& inpoel = d->Inpoel();
1460 [ + - ]: 214 : std::vector< std::vector< tk::real > > hist( d->Hist().size() );
1461 : 214 : std::size_t j = 0;
1462 [ + + ]: 282 : for (const auto& p : d->Hist()) {
1463 : 68 : auto e = p.get< tag::elem >(); // host element id
1464 : 68 : const auto& n = p.get< tag::fn >(); // shapefunctions evaluated at point
1465 [ + - ]: 68 : hist[j].resize( ncomp+1, 0.0 );
1466 [ + + ]: 340 : for (std::size_t i=0; i<4; ++i) {
1467 [ + - ]: 272 : const auto u = m_u[ inpoel[e*4+i] ];
1468 : 272 : hist[j][0] += n[i] * u[0];
1469 : 272 : hist[j][1] += n[i] * u[1]/u[0];
1470 : 272 : hist[j][2] += n[i] * u[2]/u[0];
1471 : 272 : hist[j][3] += n[i] * u[3]/u[0];
1472 : 272 : hist[j][4] += n[i] * u[4]/u[0];
1473 : 272 : auto ei = u[4]/u[0] - 0.5*(u[1]*u[1] + u[2]*u[2] + u[3]*u[3])/u[0]/u[0];
1474 : 272 : hist[j][5] += n[i] * eos::pressure( u[0]*ei );
1475 [ - + ]: 272 : for (std::size_t c=5; c<ncomp; ++c) hist[j][c+1] += n[i] * u[c];
1476 : 272 : }
1477 : 68 : ++j;
1478 : : }
1479 [ + - ]: 214 : d->history( std::move(hist) );
1480 : 214 : }
1481 : :
1482 : : // Field data
1483 [ + + ][ + + ]: 9006 : if (d->fielditer() or d->fieldtime() or d->fieldrange() or m_finished) {
[ + + ][ + + ]
[ + + ]
1484 [ + - ][ + - ]: 943 : writeFields( CkCallback(CkIndex_KozCG::integrals(), thisProxy[thisIndex]) );
[ + - ][ + - ]
1485 : : } else {
1486 : 8063 : integrals();
1487 : : }
1488 : 9006 : }
1489 : :
1490 : : void
1491 : 9006 : KozCG::integrals()
1492 : : // *****************************************************************************
1493 : : // Compute integral quantities for output
1494 : : // *****************************************************************************
1495 : : {
1496 : 9006 : auto d = Disc();
1497 : :
1498 [ + - ][ + - ]: 9006 : if (d->integiter() or d->integtime() or d->integrange()) {
[ - + ][ - + ]
1499 : :
1500 : : using namespace integrals;
1501 [ - - ]: 0 : std::vector< std::map< int, tk::real > > ints( NUMINT );
1502 : : // Prepend integral vector with metadata on the current time step:
1503 : : // current iteration count, current physical time, time step size
1504 [ - - ]: 0 : ints[ ITER ][ 0 ] = static_cast< tk::real >( d->It() );
1505 [ - - ]: 0 : ints[ TIME ][ 0 ] = d->T();
1506 [ - - ]: 0 : ints[ DT ][ 0 ] = d->Dt();
1507 : : // Compute mass flow rate for surfaces requested
1508 [ - - ]: 0 : for (const auto& [s,sint] : m_surfint) {
1509 : : // cppcheck-suppress unreadVariable
1510 [ - - ]: 0 : auto& mfr = ints[ MASS_FLOW_RATE ][ s ];
1511 : 0 : const auto& nodes = sint.first;
1512 : 0 : const auto& ndA = sint.second;
1513 [ - - ]: 0 : for (std::size_t i=0; i<nodes.size(); ++i) {
1514 : 0 : auto p = nodes[i];
1515 [ - - ]: 0 : mfr += ndA[i*3+0] * m_u(p,1)
1516 [ - - ]: 0 : + ndA[i*3+1] * m_u(p,2)
1517 [ - - ]: 0 : + ndA[i*3+2] * m_u(p,3);
1518 : : }
1519 : : }
1520 [ - - ]: 0 : auto stream = serialize( d->MeshId(), ints );
1521 [ - - ]: 0 : d->contribute( stream.first, stream.second.get(), IntegralsMerger,
1522 [ - - ][ - - ]: 0 : CkCallback(CkIndex_Transporter::integrals(nullptr), d->Tr()) );
1523 : :
1524 : 0 : } else {
1525 : :
1526 : 9006 : step();
1527 : :
1528 : : }
1529 : 9006 : }
1530 : :
1531 : : void
1532 : 8547 : KozCG::evalLB( int nrestart )
1533 : : // *****************************************************************************
1534 : : // Evaluate whether to do load balancing
1535 : : //! \param[in] nrestart Number of times restarted
1536 : : // *****************************************************************************
1537 : : {
1538 : 8547 : auto d = Disc();
1539 : :
1540 : : // Detect if just returned from a checkpoint and if so, zero timers and
1541 : : // finished flag
1542 [ + + ]: 8547 : if (d->restarted( nrestart )) m_finished = 0;
1543 : :
1544 : : // Load balancing if user frequency is reached or after the second time-step
1545 [ + + ]: 8547 : if (d->lb()) {
1546 : :
1547 : 5549 : AtSync();
1548 [ - + ]: 5549 : if (g_cfg.get< tag::nonblocking >()) dt();
1549 : :
1550 : : } else {
1551 : :
1552 : 2998 : dt();
1553 : :
1554 : : }
1555 : 8547 : }
1556 : :
1557 : : void
1558 : 8542 : KozCG::evalRestart()
1559 : : // *****************************************************************************
1560 : : // Evaluate whether to save checkpoint/restart
1561 : : // *****************************************************************************
1562 : : {
1563 : 8542 : auto d = Disc();
1564 : :
1565 : 8542 : const auto rsfreq = g_cfg.get< tag::rsfreq >();
1566 : 8542 : const auto benchmark = g_cfg.get< tag::benchmark >();
1567 : :
1568 [ + + ][ - + ]: 8542 : if ( !benchmark && (d->It()) % rsfreq == 0 ) {
[ - + ]
1569 : :
1570 [ - - ]: 0 : std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
1571 [ - - ]: 0 : contribute( meshdata, CkReduction::nop,
1572 [ - - ][ - - ]: 0 : CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
1573 : :
1574 : 0 : } else {
1575 : :
1576 : 8542 : evalLB( /* nrestart = */ -1 );
1577 : :
1578 : : }
1579 : 8542 : }
1580 : :
1581 : : void
1582 : 9006 : KozCG::step()
1583 : : // *****************************************************************************
1584 : : // Evaluate whether to continue with next time step
1585 : : // *****************************************************************************
1586 : : {
1587 : 9006 : auto d = Disc();
1588 : :
1589 : : // Output one-liner status report to screen
1590 [ + + ]: 9006 : if (thisIndex == 0) d->status();
1591 : :
1592 [ + + ]: 9006 : if (not m_finished) {
1593 : :
1594 : 8542 : evalRestart();
1595 : :
1596 : : } else {
1597 : :
1598 : 464 : auto meshid = d->MeshId();
1599 [ + - ]: 464 : d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
1600 [ + - ][ + - ]: 928 : CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
1601 : :
1602 : : }
1603 : 9006 : }
1604 : :
1605 : : #include "NoWarning/kozcg.def.h"
|
