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